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CHAPTER  I 


INTRODUCTION 
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1.1  Hydrodynamic  Stability  of  a Gas/Liquid  Interface 

Hydrodynamic  stability  is  one  of  the  important  branches  of  fluid 
mechanics.  The  original  interest  in  this  field  was  mainly  in  the  area 
of  transition  of  laminar  to  turbulent  motion.  In  recent  years,  however, 
methods  of  hydrodynamic  stability  have  been  extended  to  treat  the 
stability  of  the  interface  between  two  fluids.  The  motivation  for 
studying  dynamics  of  the  interface  comes  from  a variety  of  phenomena 
occurring  in  nature  and  in  modern  engineering  problems.  An  example  of 
the  former  which  has  fascinated  the  scientific  mind  for  a long  time  is 
the  generation  of  waves  on  the  ocean  surface  by  wind.  Since  che  wind 
motion  is  practically  incompressible,  many  of  the  earlier  investigations 
were  concerned  with  two  incompressible  fluids  in  parallel  motion.  With 
the  advent  of  space  age  it  was  felt  necessary  to  provide  some  means  of 
cooling  the  spacecraft  re-entering  the  atmosphere  at  hyper  velocities. 

An  efficient  way  of  achieving  this  is  to  inject  a liquid  coolant  through 
the  nose  of  the  spacecraft  so  that  a layer  of  liquid  is  maintained  on 
the  windward  surface.  Thus,  effective  cooling  can  be  accomplished  if 
the  liquid  film  stays  adhered  to  the  surface.  Therefore  stability  of 
liquid  films  in  high  speed  gas  environments  received  a great  deal  of 
attention  in  the  1970's. 

It  is  a fact  of  experience  that  wind  blowing  over  water  generates 
waves  on  the  interface  and  that,  under  certain  conditions,  these  waves 
grow  in  size.  Once  a wave  reaches  a significant  height  droplets  are 
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stripped  off  from  the  wave  crests  and  get  entrained  in  the  air.  Hence 
the  liquid  is  lost  not  only  due  to  the  simple  process  of  evaporation 
but  also  due  to  the  entrainment.  In  fact,  current  investigations  show 
that  the  latter  mechanism  is  more  dominant. 

It  is  the  purpose  of  the  present  work  to  introduce  the  element  of 
interface  mass  transfer  into  the  stability  problem.  Several  simplifying 
assumptions  have  been  made  to  gain  a first  insight  into  the  nature  of 
the  problem  and  to  permit  analytical  tractability . To  this  end  only 
the  linear  stability  problem  is  considered.  Physically,  this  assumption 
implies  that  the  waves  on  the  interface  are  of  infinitesimal  amplitude 
or,  put  differently,  the  wave  amplitude  is  much  smaller  than  the  wave- 
length. It  is  unlikely  that  the  mechanism  of  entrainment  would  be 
significant  under  these  circumstances  and  hence  only  the  mass  transfer 
due  to  evaporation  is  taken  into  account.  Another  key  assumption 
introduced  in  the  present  study  is  that  the  gas  motion  is  incompressible. 
It  may  be  argued  that  mass  transfer  effects  would  be  negligible  in  the 
incompressible  case  and  therefore  it  is  the  case  of  compressible  gas  that 
is  worth  examining.  However,  in  many  natural  processes  and  engineering 
problems  the  air  flew  is  incompressible;  for  example,  interactions 
between  wind  and  ocean  surface,  dispersion  of  pollutants  into  the 
atmosphere  from  large  bodies  of  water,  two-phase  flows  and  industrial 
drying  processes.  Apart  from  these  practical  applications,  the  fact 
that  the  incompressible  problem  has  not  yet  been  completely  understood 
provides  motivation  for  the  present  work.  Also,  the  experience  gained 


in  Che  study  of  boundary  layer  stability  has  shewn  that  the  compressible 
problem  is  considerably  more  complicated  than  its  incompressible  counter^ 
part.  Indeed,  a complete  spectrum  of  eigenvalues  of  the  simple  Blasius 
flat  plate  boundary  layer  was  obtained  only  very  recently.  Thus  an 
examination  of  the  incompressible  case  appears  to  be  the  logical  first 
step.  Under  the  assumption  of  incompressible  gas  and  the  restriction 
of  evaporative  mass  transfer  the  rate  of  mass  injection  into  the  gas  is 
expected  to  be  small.  It  will  be  shown  that  under  these  conditions  the 
exponential  steady-state  velocity,  temperature  and  concentration  profiles 
reduce  to  linear  profiles.  In  addition  to  the  assumptions  outlined  above 
the  classical  assumption  of  parallel  mean  (or  steady- state)  flew  is  also 
made.  It  is  then  possible  to  solve  the  governing  equations  in  a closed 
form. 

The  present  work  is  an  attempt  to  determine  a qualitative  and 
quantitative  estimate  of  the  influences  of  mass  transfer  on  the  liquid 
film  stability  problem.  This  endeavor  is  pursued  with  the  intention  of 
producing  a useful  analytical  framework  whi>  Li  can  treat  such  interesting 
aspects  as  (i)  amplification  and  phase  velocity  curves  (ii)  neutral 
stabili.y  curves  (iii)  stress,  temperature  and  concentration  pertur- 
bations at  the  interface  and  (iv)  energy  transfer  mechanisms.  The 
analysis  presented  here  departs  from  the  customary  assumption  of 
neglecting  the  instabilities  in  gas,  which  amounts  to  ignoring  the 
eigenvalue  in  gas  disturbance  equations.  The  consequences  of 
relaxing  this  assumption  are  carefully  examined. 
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1.2  Review  of  Literature 

It  has  been  recognized  for  a long  time  that  waves  are  Initiated 
on  a liquid  surface  due  to  some  kind  of  instability.  The  growth  or 
decay  of  these  waves  depends  on  whether  energy  is  transferred  to  or 
removed  from  them.  Hence  the  study  of  mechanisms  of  energy  transfer 
has  received  considerable  attention  in  the  literature.  Ursell^  pro- 
vides an  excellent  survey  of  wind  generated  waves  on  deep  water.  A 

2 3 

widely  respected  theory  of  wave  initiation  is  due  to  Miles  ’ . He 
proved  that  the  rate  of  energy  transfer  to  a wave  of  speed  c is  pro- 
portional to  the  profile  curvature  - U"(c)  in  the  gas.  This  theory 

predicts  exponential  wave  growth  whereas  experiential  observations 

4 

show  that  the  initial  wave  growth  is  linear.  It  was  Phillips  who 
explained  the  initial  linear  wave  growth  by  proposing  a resonance 
mechanism  between  turbulent  pressure  fluctuations  and  interface  dis- 
turbances. A combined  Miles-Phillips^  theory  provides  a very  good 
qualitative  description  of  wave  initiation  and  growth.  Lighthill^ 
has  offered  a remarkable  physical  interpretation  of  Miles'  theory. 

There  are  three  principal  types  of  instabilities  which  are  of 
interest  in  interface  stability  problems.  These  are  discussed  below 
and  the  relevant  literature  is  reviewed. 

(i)  Tollmien-Schlichting  Instability: 

The  Tollmien-Schlichting  mechanism  has  been  studied  extensively 
in  connection  with  the  stability  of  laminar  boundary  layers  (Ref.  7). 


r 


6 


This  instability  is  the  result  of  continual  energy  transfer  from  the 

8 9 10  11 

mean  flow  to  the  disturbance.  Benjamin  ’ , Landahl  and  Skripachev 
investigated  the  stability  of  laminar  boundary  layers  over  flexible 
surfaces.  Their  analyses  predicted  three  different  stability  modes, 
viz.  modified  forms  of  the  Tollmien-Schlichting  mode,  gravity  waves 


(body  force  instabilities)  and  Kelvin-Helmholtz  instability.  The 

latter  two  types  will  be  described  shortly.  In  recent  years  the  higher 

stability  modes  associated  with  the  incompressible,  laminar  flat  plate 

12  13 

boundary  layer  have  been  obtained  by  Gastor  and  Jordinson  and  Mack 

Mack's  calculations  show  that  the  number  of  eigenvalues  is  finite  for 

a velocity  profile  which  is  sufficiently  smooth  at  the  outer  edge  of 

the  boundary  layer.  This  work,  along  with  the  analysis  of  DiPrima  and 

Habetler^,  shows  that  in  a finite  interval  there  exists  an  infinite 

discrete  eigenvalue  spectrum.  Gallagher  and  Mercer and  Deardorff 

examined  the  eigenvalue  spectrum  of  a plane  Couette  flow  and  concluded 

that  it  is  stable  to  infinitesimal  disturbances.  The  stability  of  a 

plane  Couette-Poiseuille  flow  with  uniform  cross  flow  has  been  inves- 
18 

tigated  by  Haines 

(ii)  Instability  due  tc  shear  and  pressure  perturbations  at  the  inter- 
face: 

Small  wavy  disturbances  on  the  interface  become  unstable  when 
pressure  and  shear  perturbations  overcome  the  stabilizing  effect  of 
gravity  and  surface  tension.  A special  case  of  this  instability 
mechanism  is  the  well-known  Kelvin-Helmholtz  instability  (Ref.  19) 


between  two  incompressible,  inviscid  fluids  in  parallel  motion.  Miles 


7 


generalized  this  simple  case  to  parallel  shear  flows.  In  a classic 
21 

paper  Benjamin  presented  results  for  pressure  and  shear  stress  per- 
turbations exerted  by  a viscous  incompressible  fluid  on  a wavy  wall. 

He  considered  the  cases  of  rigid,  flexible  yet  solid,  and  a completely 
mobile  wavy  wall.  In  the  course  of  his  investigation  he  also  laid 

dcwn  the  requirements  under  which  a moving  wavy  interface  could  be 

22 

approximated  as  rigid.  These  results  were  applied  by  Craik  to 

wind-generated  waves  in  thin  liquid  films.  He  discovered  that  the 

dominant  instability  mechanism  in  thin  films  is  due  to  shear  pertur- 
23 

bations.  Lighthill  analyzed  the  shear  and  pressure  perturbations 

exerted  by  a viscous  compressible  gas  on  a rigid  wavy  wall  neglecting 

the  effects  of  heat  and  mass  transfer  at  the  wall.  This  work  was 
24 

amplified  by  Inger  to  include  heat  and  mass  transfer  effects. 

(iii)  Rayleigh-Taylor  instability: 

This  instability  is  due  to  body  forces  pointing  from  the  heavier 

fluid  to  lighter  fluid  (Ref.  25).  The  problem  of  liquid  film  stability 

26 

in  a body  force  field  was  solved  by  Nayfeh  and  Saric 

Several  investigators  have  studied  the  linear  stability  of 
liquid- gas  interface  under  various  conditions  and  a concise  summary  of 
their  works  is  provided  in  Table  I at  the  end  of  this  chapter.  It 
should  be  noted  that  the  liquid  motion  is  treated  as  incompressible  and 
laminar  in  all  cases,  and  effects  of  mass  transfer  are  neglected  unless 
mentioned  otherwise.  Major  features  of  the  present  work  are  also  listed 


in  this  table  to  facilitate  comparison  with  other  investigations. 
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1.3  Outline  of  Present  Work 

It  was  pointed  out  in  the  previous  section  that  steady-state  expo- 
nential mass  transfer  profiles  reduce  to  linear  profiles  without  mass 
transfer  when  the  evaporation  rate  is  small.  This  suggests  that  the 
zero  mass  transfer  problem  with  linear  profiles  be  solved  first.  Since 
the  zero  mass  transfer  problem  is  relatively  simpler,  the  experience 
gained  in  its  formulation  would  be  very  useful.  Also,  this  problem  can 
be  used  to  evaluate  the  importance  of  the  often-made  assumption  of 
neglecting  gas  instability.  The  formulation  of  the  zero  mass  transfer 
problem  is  described  in  Chapter  II  and  its  solution  is  presented  in 
Chapter  IV.  Two  methods  of  solution  have  been  used  - (i)  a small  pertur- 
bation scheme  which  yields  one  eigenvalue  in  the  long  wavelength  approx- 
imation and  (ii)  exact  analytical  solution  for  arbitrary  wavelengths. 

The  mass  transfer  formulation  is  considered  in  Chapter  III  and 
follows  closely  the  procedures  of  Chapter  II.  Chapter  V is  concerned 
with  the  solution  of  the  mass  transfer  problem.  All  the  important  steps 
of  mathematical  analysis  are  incxuded  in  Chapters  II  through  V and  the 
details  are  confined  to  the  Appendices.  This  was  done  in  order  to 
preserve  clarity  of  the  presentation.  The  results  of  numerical  com- 
putations are  given  in  the  form  of  amplification  and  phase  velocity 
curves  in  Chapter  VI.  The  conclusions  derived  from  the  present  investi- 
gation are  summarized  in  Chapter  VII. 


surface  for  the  onset  of  instability. 
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the  wave  slope  and  pressure  perturbation 
in  phase  with  the  wave  height. 


elude  mass  transfer  effects 


Kotake  39  Incompressible,  Parallel  flow.  Analyzed  the  steady-state  problem  and 

(1973)  laminar,  infinite  depth  found  that  interface  evaporation  leads 
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CHAPTER  II 


FORMULATION  OF  THE  ZERO  MASS  TRANSFER  PROBLEM 


17 


2.1  The  Purpose  of  Solving  the  Zero  Mass  Transfer  Problem 

1.  It  was  mentioned  in  Sec.  1.2  that  different  types  of  mecha- 
nisms have  been  proposed  to  explain  the  transfer  of  energy  to  the 
surface  waves.  These  mechanisms  are  briefly  summarized  below. 

(i)  Energy  transfer  from  the  mean  velocity  profile  in  the  gas: 

2 21 

This  mechanism  was  proposed  by  Miles  and  Benjamin  . The 
critical  point  (y-location  at  which  the  disturbance  phase  speed 
equals  the  u-velocity  component)  of  the  velocity  profile  lies 
within  the  gas  boundary  layer.  Energy  is  fed  from  the  mean  flow 
to  the  interface  disturbance  continually  with  time  resulting  in  a 
Tollmien-Schlichting  type  of  instability.  However,  the  experi- 
ments of  Cohen  and  Hanratty^  and  Plate,  et  al.^  indicate  other- 
wise. 

(ii)  Energy  transfer  from  the  mean  velocity  profile  in  the  liquid: 
The  critical  point  lies  inside  the  liquid  layer  and  again  the 

energy  transfer  to  the  interface  occurs  as  a result  of  Tollmien- 

Schlichting  instability.  This  mechanism  was  investigated  by  Feld- 
29 

man  and  Miles  for  large  liquid  Reynolds  numbers.  In  this  mode 

the  phase  speed  c^.  is  less  than  the  interface  velocity  (c^  < u^j.) 

and  these  waves  are  sometimes  referred  to  as  'slow  waves.'  These 

22 

waves  have  been  observed  in  the  experiments  of  Craik  and  Saric 
3 A 

and  Marshall J . However,  Cohen  and  Hanratty  reject  this  mode  of 

energy  transfer  because  they  did  not  observe  sloT»  waves  in  their 

experiments.  One  reason  may  be  that  they  did  not  use  sufficiently 

34 

small  thicknesses  in  their  work.  Saric  and  Marshall  mention 


t.  ’t  the  occurrence  of  slow  waves  may  be  due  to  nonlinear  effects. 

(iii)  Energy  transfer  due  to  pressure  and  shear  perturbation 

exerted  by  gas  on  the  disturbed  interface: 

The  classical  Kelv in-Helmholtz  stability  is  a special  case 

30 

of  this  mechanism.  Cohen  and  Hanratty  proposed  that  this  is  the 

sole  mechanism  of  energy  ‘..ansfer.  They  found  that  for  'fast' 

waves  (critical  point  inside  the  gas  or  > u^)  or  'thick'  films 

the  component  of  pressure  perturbation  in  phase  with  the  wave  ele- 

22 

vation  are  dominant.  Craik  , however,  discovered  that  for  'slow' 
waves  on  very  thin  films,  the  pressure  perturbation  component  in 
phase  with  the  wave  elevation  and  shear  perturbation  component  in 
phase  with  the  wave  slope,  are  dominant. 

Another  model  based  on  energy  transfer  due  to  pressure  vari- 
ations at  the  interface  is  Jeffrey's^  sheltering  hypothesis.  It  is 
improoable  that  the  mechanism  of  sheltering  (i.e.  drag  force  exerted 
on  a wave  due  to  flow  separation  near  the  wave  crest)  will  be  impor- 
tant in  the  case  of  waves  with  amplitudes  very  small  compared  to 
their  wavelengths.  The  latter  is  assumed  in  all  the  linear  analyses 
including  the  present  work. 

(iv)  Energy  transfer  due  to  a resonance  mechanism  between  turbulent 

pressure  fluctuations  and  surface  disturbances: 

4 

Proposed  by  Phillips  , this  mechanism  is  believed  to  be  respon- 
sible for  initiation  of  short  waves  on  the  interface.  Cohen  and 
30 

Hanratty  rule  out  this  mode  on  the  basis  that  a smooth  liquid  sur- 
face is  'observed'  even  in  the  presence  of  turbulent  flow  (not  a very 


19 


31 

convincing  argument  at  all!).  Plate,  et  al.  obtained  turbulent 
pressure  fluctuations  at  the  liquid  surface  indirectly  through  the 
measurement  of  longitudinal  velocity  fluctuation  u'.  Their  experi- 
ments come  closest  to  verifying  Phillips'  theory.  They  conclude 
that  the  resonance  mechanism  does  not  appear  to  be  significant  and 
this  may  be  a justification  for  neglecting  turbulence  interactions 
in  the  present  work.  One  of  the  aims  of  solving  the  zero  mass  trans- 
fer problem  is  to  try  to  shed  some  light  on  the  mechanisms  of  energy 
transfer. 

2.  As  seen  earlier  in  Sec.  1.2  some  investigators  have  a priori 

33 

neglected  the  instability  in  the  gas.  For  instance,  Nachtscheim 
35 

and  Starkenberg  observe  thac  c/ug<<  1 for  inviscid  supersonic 

external  flow,  where  c is  the  phase  velocity  of  the  disturbance 

37 

and  ug  is  the  gas  velocity.  Bordner  presents  an  order  of  magni- 
tude argument  for  neglecting  c in  the  gas  disturbance  equations  when 

the  external  flow  is  viscous  and  compressible.  In  his  work  on  cross- 
24 

hatching  Inger  also  assumes  that  the  interface  behaves  like  a rigid 


22 

wavy  wall  relative  to  the  gas  flow.  Craik  analysed  thin  liquid 


films  using  the  pressure  and  shear  perturbations  derived  by  Benjamin 
for  a rigid  wavy  boundary.  This  amounts  to  having  the  critical 
point  located  at  the  wavy  boundary  and  it  also  means  that  the  phase 
speed  c is  negligible  in  the  gas  disturbance  equations.  It  ought  to 
be  mentioned  at  this  point  that  the  assumptions  (a)  the  phase  speed 
relative  to  gas  speed  is  negligible  (b)  the  interface  is  steady  and 
rigid  for  gas  disturbance  equations  and  (c)  the  critical  point  in  the 


21 


gas  is  located  at  the  wavy  boundary,  are  all  equivalent. 

22 

Craik  indicates  that  even  for  a viscous  incompressible  gas 
c/ue<<l  is  sufficient  to  make  the  above  assumption  (which  is  under- 
standable in  the  case  of  an  inviscid  compressible  gas).  It  is  ques- 
tionable how  Craik' s assumption  (i.e.  c/ue<<l)  satisfies  the  require- 

21 

ment  laid  down  by  Benjamin  that 


me 

u'(0)  <C  1 


(2.1.1) 


in  order  to  make  the  rigid  wavy  wall  assumption.  In  the  last  inequa- 
lity c = c/ug,  u' (0)  is  the  non-dimensional  slope  of  the  velocity 
profile  evaluated  at  wall  and 


m = [aRu' (o) ] 


(2.1.2) 


where  R is  the  gas  Reynolds  number  based  on  some  characteristic 
thickness  and  a = 2ir/A,  A being  the  wavelength  of  the  rigid  wavy  wall. 
Thus  Benjamin's  criterion  requires  that 


[oRu'(0)]1/3  «1 

u ' (0 ) 


(2.1.3) 


It  is  clear  from  Eq.  (2.1.3)  that  for  a low  speed  (incompressible) 
gas  with  moderate  Reynolds  number  and  sufficiently  large  a(small 
ripples  on  the  interface)  Benjamin's  criterion  may  not  hold.  Hence 
c/ug  may  be  very  small  compared  to  unity  but  still  Eq . (2.1.1)  could  be 
violated  and  consequently,  a rigid  wavy  interface  assumption  cannot 
be  made.  An  illustrative  example  will  be  considered  in  Chapter  VI. 


21 


In  the  present  analysis  the  external  gas  is  viscous  and  incom- 
pressible and  hence  the  rigid/steady  interface  assumption  is  not  made. 
Hence  phase  speed  terms  appear  in  both  the  gas  and  liquid  disturbance 
equations.  The  stability  problem  is  solved  both  with  and  without  this 
assumption,  thus  providing  a method  of  checking  its  validity. 

3.  An  incompressible,  viscous,  laminar  flow  of  gas  with  constant 
properties  over  an  incompressible  laminar  viscous  liquid  is  considered. 
A turbulent  velocity  profile  in  the  gas,  however,  can  be  treated  as  a 
'quasi-laminar ' profile  with  augmented  viscosity.  The  assumption  of 
laminar  flow  simplifies  the  analysis  greatly  but  it  is  not  very 
realistic.  In  the  turbulent  case  it  would  imply  that  the  laminar 
sublayer  be  large  compared  to  the  disturbance  wavelength  — a require- 
ment rarely  met  in  practice. 

The  zero  mass  transfer  problem  is  solved  for  a linear  steady 


state  velocity  profile  in  both  gas  and  liquid.  The  linear  profile  in 
the  gas  can  be  justified  to  some  extent  in  the  laminar  case  but  it  may 


be  a poor  approximation  in  the  turbulent  case.  The  laminar  flow  and 
linear  velocity  profile  assumption  can  be  justified  for  the  liquid 
since  the  liquid  Reynolds  number  is  usually  small  and  since  Craik's 
experiment  confirms  a linear  profile.  The  linear  velocity  problem  is 
solved  in  the  present  work  due  to  the  following  reasons: 

(i)  the  Orr-Sommerf eld  equation  has  an  exact  solution. 


(ii)  the  linear  velocity  profile  is  the  simplest  viscous  profile 
which  represents  a physically  possible  flew. 

(iii)  the  linear  stability  of  plane  Couette  flow  has  been  extensively 


22 


studied  (e.g.  Gallagher  and  Mercer"^,  Deardorf f ^ and  it  has  been 
found  to  be  unconditionally  stable. 

(iv)  The  exponential  steady-state  velocity  profiles  with  mass  trans- 
fer reduce  to  linear  profiles  for  small  (but  non-zero)  rates  of  mass 
transfer.  This  last  reason,  in  particular,  provided  motivation  for 
studying  the  zero  mass  transfer  problem  with  linear  steady-state 
profiles. 


2.2  The  Steady-State  Problem 


As  mentioned  in  Sec.  1.1  the  steady-state  or  the  mean  flow  is 
assumed  to  be  incompressible  and  parallel  (i.e.  — = 0 and  v^  = = 0) . 

The  liquid  motion  is  assumed  to  be  a plane  Couette  flow  and  hence  it  is 
an  exact  solution  of  the  full  Navier-Stokes  equations.  The  gas  motion, 
on  the  other  hand,  is  assumed  to  be  a boundary  layer  flow  which  is 
approximately  parallel. 

Let  h and  6 be  the  height  of  the  liquid  Couette  flow  and  the 
boundary  layer  thickness  respectively  (Fig.  la).  The  gas  Prandtl 
number  P^  is  assumed  to  be  unity  so  that  the  thermal  and  velocity 
boundary  layer  thicknesses  are  identical . Let  subscripts  1 and  2 
denote  the  liquid  and  gas  respectively.  The  steady-state  governing 
equations  are  — 

Liquid : 


1.  x-momentum 


r 


23 


| 


2.  y-momentum 


dPl 

dT=  ~ pi8 


3.  Energy 


Gas ! 

4 . x-momentum 

5 . y-momentum 


6.  Energy 


= 0 


n 


2 


0 


(2.2.3) 


(2.2.4) 


(2.2.5) 


(2.2.6) 


These  equations  show  that  all  the  flow  quantities  vary  only  with 
respect  to  the  vertical  co-ordinate  y.  In  Eq.  (2.2.6)  the  viscous 
dissipation  term  is  neglected  consistent  with  the  assumption  of 
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incompressibility.  It  may  also  be  mentioned  here  that  continuity 
equations  for  the  gas  and  liquid  are  identically  satisfied  due  to  the 
parallel  flow  assumption.  Eqs.  (2.2.1)  - (2.2.6)  are  six  equations 
in  six  unknowns  uj , U2,  pi , P2>  T]»  T2 • The  combined  order  of  this 
system  is  10  and  an  equal  number  of  boundary  conditions  is  required 
for  a unique  solution.  These  conditions  are  — 

1.  No  slip  at  the  wall 

uj(-h)  = 0 (2.2.7) 


2.  Boundary  layer  edge  condition  on  velocity 


u2(6)  = u 


(2.2.H) 


3.  No  slip  in  tangential  velocities  at  the  interface 


u 1 (0)  = u2 


(2.2.9) 


4.  Balance  of  shear  stresses  at  the  interface 


iui 

dy 


dun 


y = 0 ^2dy 


y = 0 


(2.2.10) 


5.  Constant  temperature  or  adiabatic  wall 


Ti(-h)  = Tw  or 


dT  1 
dy 


= 0 


y=-h 


(2.2.11) 


6.  Boundary  layer  edge  condition  on  temperature 


T2(6)  = Te 


(2.2.12) 


A 


7.  Energy  balance  at  the  interface 

dT,  I 


dT9 

y = 0 dy  y = 0 

(2.2.13) 

8.  No  jump  in  temperature  at  the  interface 

Ti<0)  = T2(o> 

9.  Boundary  layer  edge  condition  on  pressure 

P2(6)  = Pe 

10.  Balance  of  normal  stresses  at  the  interface 

Pl(0)  = P2(0) 


(2.2.14) 


(2.2.15) 


(2.2.16) 


The  solution  of  Eqs.  (2.2.1)  - (2.2.6)  subject  to  (2.2.7)  - 
(2.2.16)  is 

Liquid  velocity  profile 


Gas  velocity  profile 


Liquid  pressure  profile 


!£  i 1 * h 
Ue'  6 , ^2  h 

+ U1  6 


- H h + i 

u p 6 6 

_2  = 

Ue 

Ul  <-> 


p 1 = p - Pi  gy 


(2.2.17) 


(2.2.18) 


(2.2.19) 
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Gas  pressure  profile 


Liquid  temperature  profile 


(2.2.20) 


1 - 


_1  h 

k^  6 


T T k. 

w w _2_  ri 

T T k.  6 

g_...  1 + _£ L_ 

k„  . h k„ 


constant  temperature  wall 


i4i 

kx  6 


1 + J2*L 

kl  6 


= 1 adiabatic  wall 


(2.2.21) 


(2.2.22) 


Gas  temperature  profile 


1 - 


w 


1 +11  Jl 

1 + k1  6 


T k.  , 

w £ 

Te  kl  5 

1 h 
1 kx  6 


constant  temperature  wall  (2.2.23) 


= 1 adiabatic  wall 

Interface  quantities  are  obtained  from  Eqs.  (2.2.17)  - (2.2.24) 
putting  y = 0.  Thus 
Interface  velocity 


(2.2.24) 

by 


if 


^2h 

hi  « 


(2.2.23; 


1 + — 
vl 


2 h 


Interface  pressure 


if 


(2.2.26) 


The  Unsteady  Problem  — Governing  Equations 


When  the  steady  interface  configuration  of  Sec.  2.2  is  disturbed 
the  resulting  unsteady,  two  dimensional,  incompressible  motion  is 
governed  by  the  following  equations: 

Liquid : 

1.  Continuity 

3~u,  Sv. 

-A  + -A  = 0 (2.3.1) 
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4.  Energy  (static  enthalpy  form) 


3h^  3h.^  3p^  3p^  3p 

pi  l^r  + “i  + \ 3tJ  = ~ + “i  37~ + vi  37- 


'32T.  32T  ' 

+ k.  + i + p. 


^(3x2  3y2 


(2.3.4)’ 


Combining  this  equation  with  the  equation  of  state  for  the  liquid 

hl  hl  Tp 

and  neglecting  viscous  dissipation  and  pressure  gradient  terms  it  can 
be  shown  that  (Appendix  A) 


3T1  _ 3T  3T  [32T  32T 

PiC  , T—  + U,  TT-  + V,  T—  = k,  + 


i pi  i^t  t“i  ax  - vi  wr  ij^rT^r 


(2.3.4) 


5.  Continuity 


3u„  3v„ 

_£  + _£  = ° 
3x  3y 


(2.3.5) 


6.  x-momentum 


3U2,,  k2.  1 8p2  . ..  [3S  . 8S 


at  2 ax  2 ay  p9  ax  v2  2 

£ aX 


(2.3.6) 
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7.  y- momentum 


3v2  _ 3v2  — 3v2 

at-  + u2  "ST  + V2  ■57" 


i_!!i  + v f!S  + 

P2  8y  "2  8x2  + 8y2 


8.  Energy  (static  enthalpy  form) 


^2  _ ^^2  — 3^2  3p2  — 3p2  ^P2 

P2  JtT  + U2  87”  + V2  87",  = 3t“  + U2  8x~  + V2  TT 


82T2  82T2' 

+ k„  + + y,4>_ 

2 8x2  8y2  22 


combining  this  equation  with  the  relation  for  an  ideal  gas 


(8h/8T)  = c 

P P 


and  neglecting  pressure  gradient  and  viscous  dissipation  terms  the 
result  is 


'8T  8T  8T  ] f82T  82T 

P~c  9 + u~  t — + v.  - — = k„  + 

2 P2|8t  2 8x  2 8y  ^ 2^  3y 2 


It  is  assumed  that  the  unsteady  motion  is  confined  to  the  region 
-h  y < 5,  i.e.  the  steady-state  gas  boundary  layer  thickness  6 is 
unaltered.  This  is  justifiable  in  view  of  the  assumption  of  small 


perturbation  motion  to  be  introduced  later  (Sec.  2.5). 
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2.4  The  Unsteady  Problem  — Boundary  Conditions 

Eqs.  (2.3.1)  - (2.3.8)  form  a system  of  elliptic  partial 
differential  equations  requiring  boundary  conditions  to  be  specified 
along  the  entire  boundary  of  the  domain  in  (x,y)  plane.  However,  it 
will  be  apparent  in  Sec.  2.6  that  the  unsteady  perturbations  are  sinus- 
oidal with  respect  to  the  x co-ordinate  and  hence  bounded.  Thus  boun- 
dary conditions  need  be  specified  only  along  the  y co-ordinate  direction 
and  along  the  interface.  In  this  section  the  appropriate  boundary 
and  interface  conditions  are  developed. 

1.  No  slip  at  the  wall 

u1(x,y,t)  =0  , y = -h  (2.4.1) 

2.  Boundary  layer  edge  condition  on  the  u-velocity  component 


u2(x,y,t)  = ug  , y = 6 


(2.4.2) 


3.  No  slip  in  the  tangential  velocity  at  the  interface  (Fig.  2) 

Let  VD  and  V be  the  liquid  and  gas  velocity  vectors  at  point 
1 A2 

p relative  to  the  interface.  The  condition  of  no  slip  requires  that 

the  components  of  V and  vT  along  the  interface  be  continuous,  i.e. 
R1  R2 


ds 


• ds 


where  ds  is  a directed  line  segment  of  the  interface  at  point  p.  Now, 
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and 


- V 


if 


where  and  are  liquid  and  gas  velocity  vectors  respectively 
relative  to  a stationary  observer.  V ^ is  the  interface  velocity 
vector  with  respect  to  a stationary  observer.  Hence, 


or 


Since 


i.e. 


(V1  “ Vif}  ' ds  = (V2  " Vif}  ' ds 


• ds  = • ds 


vi  = ^Ui,vi^  ’ V2  = ^u2,v2^  and  ds  = (dx»dy> 


u^dx  + v^dy  = u^dx  +v^dy 


(u2  - Ul)  = (vx  - v2)  g 


If  y = n(x,t)  is  the  equation  describing  the  unsteady  interface,  the 
above  relation  gives 


u2  " ui  = (vl  - V U at  y = n(x’° 


(2.4.3) 


'"■T"  “ 


~~ 
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For  a flat  interface  nx  = 0 and  (2.4.3)  reduces  to  (2.2.9). 

4.  Balance  of  shear  stresses  at  the  interface 

Ti  = t2  on  y = 

it  can  be  shown  by  considering  the  equilibrium  of  a triangular 
element  at  the  interface  (Appendix  B)  that 


^x2 

T = M 

l+n  2 


’ ' 
3u  . 3v 

2unx 

3u 

' 

3v 

1 

h 

1 

[CT5 

l+n  2 

3x 

3y; 

Thus  shear  balance  requires  that 


1-n  2 


l+n  2 

X 


M2 


3y 


1-n 


l+n 


auj 

3y 


+ !^1 
9x 

J 

2u2\  J 

( — 

3v 

l+n  2 ^ 

X 

[3x 

1 ^ 
< 1 
I*-1 

2Mlnx 

[9ul 

9x  j 

l+n  2 

>3x 

3y  , 

on  y = n(x,t) 


(2.4.4) 


For 

a flat  interface  n 

X 

= °>  \ = ^2 

= 0 and  Eq . 

(2.4.4)  reduces  to 

Eq. 

(2.2.10). 

5. 

Constant  temperature 

wall 

T^x.y.t)  = 

H 

C 

II 

1 

or 

Adiabatic  wall 

3T.  (x,y,t) 

(2.4.5) 
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6.  Boundary  layer  edge  condition  on  temperature 


T2(x,y,t)  = Te,  y = 6 


7.  Energy  balance  at  the  interface 

3T  3T2 


(2.4.6) 


or 


klVTl  ‘ n = k2VT2  * n 


where  n is  the  unit  normal  to  the  interface  at  any  point.  Now,  if  the 
interface  equation  is  written  as 

F(x,y,t)  = y - n(x,t)  = 0 


then, 


- VF 

n = | VF  r 


and  henae 


klVTl 


VF 


k^VT^  " VF  on  F = 0 or  y = n(x,t)  (2.4.7) 


Expanding  (2.4.7) 


9T1  3F  ^1  3F 

3x  3x  3y  3y 


= k„ 


^2  3F  t dT2  3F 
3x  3x  3y  3y 
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using  the  relation  F = y - o(x,t). 


3T 

3T1 

= K 

3T2 

3] 

3x  nx 

3y  . 

2 

3x 

3y  - 

For  a flat  steady  interface  H 0,  nx  = 0 
reduces  to  Eq.  (2.2.13) 

8.  No  temperature  jump  at  the  interface 


the  above  equation 


T^(x,y , t)  = T2(x,y,t)  on  y = n(x,t)  (2.4.8) 


9.  Boundary  layer  edge  condition  on  pressure 


P2(x,y,t)  = Pe  at  y = 6 (2.4.9) 

10.  Balance  of  normal  stresses  at  the  interface 

Referring  to  Fig.  3a  it  is  seen  that  the  discontinuity  in  the 
normal  stresses  across  the  interface  must  be  balanced  by  surface 
tension.  Thus 


°2  ' °1  = £ on  y 


n(x,t) 


where  r is  the  surface  tension  and  R is  the  radius  of  curvature.  Now, 
for  a surface  that  is  concave  downward,  the  curvature  is  given  by 


1 

R 


-n 


XX 


<l+n*> 


3/2 


Also,  as  shown  in  Appendix  B,  the  expression  for  normal  stress  is 
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a 


P + 


_2U_ 

l+n2 

X 


3u 

3x 


2unx 


3u  3v 
3y  3x 


Finally,  the  normal  stress  condition  becomes 


(p2  - p1) (1+n2)  - 2n2 


3u2 

3U1 

3v2 

a? 

P2  IbT  ” 

P1  3jT 

4 

- 2 

M1  3y  j 

+ 2n 


x 


at  y = n(x,t)  (2.4.10) 

In  the  steady  state  case  Eq.  (2.4.10)  reduces  to  Eq.  (2.2.16) 

The  boundary  conditions  developed  so  far,  (2.4.1)  through 
(2.4.10),  are  the  same  as  those  for  the  steady-state  problem,  viz. 
(2.2.7)  through  (2.2.16).  In  the  steady-state  case  there  were  six 
unknowns  u^,  u2,  p^,  P2>  T^»  and  T2*  In  unsteady  case,  however, 
there  are  eight  unknowns,  u^ , u2,  , "p2,  T^,  T2,  v^  and  v,,.  Thus  v^ 

and  v2  are  the  two  additional  unknowns  which  appear  as  second 
derivatives  with  respect  to  x and  y and  also  as  first  derivatives  with 
respect  to  time  (Eq.  (2.3.3)).  As  will  be  shown  in  Sec.  2.6  the  x and 
t dependence  can  be  eliminated  by  assuming  a suitable  form  of  interface 
configuration.  Thus  one  need  be  concerned  only  with  the  boundary 
conditions  on  v^  and  v2  in  the  y direction.  Since  v^  and  v2  appear  as 
second  derivatives  w.r.t.  y in  Eqs.  (2.3.3)  and  (2.3.7),  four 


conditions  are  required  on  and  v These  conditions  are,  the  wall 
condition  on  v, , the  boundary  layer  edge  condition  on  v^  and  two 
interface  matching  conditions.  These  four  conditions  are  derived  below. 

11.  No  penetration  condition  at  the  wall 

In  the  absence  of  any  mass  transfer  through  the  wall  the  no 
penetration  condition  is 

v^x.y.t)  = 0 , y = -h  (2.4.11) 

12.  Boundary  layer  edge  condition  on  the  vertical  velocity  component  v ^ 

The  flow  configuration  of  Fig.  1 behaves  as  if  it  is  bounded 
between  two  walls  at  y = -h  and  y = 6 even  in  tue  disturbed  state,  and 
this  would  require  that  the  v^  be  zero  at  y = 6,  i.e. 

v2(x,y,t)  = 0 , y = 6 (2.4.12) 

13.  Kinematic  condition  on  v^  at  the  interface 

Since  there  is  no  mass  transfer  across  the  interface,  the  no 
penetration  condition  (commonly  referred  to  as  the  kinematic  condition)  at 
the  deformable  interface  reads 

D1F 

jj— - = 0 on  F = 0 or  y = n(x,t) 


where 


(2.4.13) 


37 


14.  Kinematic  condition  on  v2  at  the  interface. 

Similar  reasoning  as  above  gives  the  no  penetration  condition  on 
the  gas  side  as 


where 


D2F 

= 0 or  F = 0 or  y = n(x,t) 


D2t  3t  + U2  3x  + V2  3y 


(2.4.14) 


In  the  steady-state  case  the  boundary  conditions  (2.4.11)  through 
(2.4.14)  are  identically  satisfied. 

Eqs.  (2.3.1)  through  (2.3.8)  supplemented  by  boundary  conditions 
(2.4.1)  through  (2.4.14)  complete  the  general  formulation  of  the 
unsteady  problem.  It  is  noted  that  this  problem  is  highly  nonlinear. 


2.5  Small  Perturbation  Formulation  of  the  Unsteady  Problem 

The  solution  of  the  unsteady  problem  in  its  most  general 
nonlinear  form  presents  a formidable  task  and  hence  a small  perturbation 
solution  is  attempted.  The  unsteady  motion  is  viewed  as  a small  per- 
turbation on  the  steady-state  problem  of  Sec.  2.2.  Accordingly,  every 
dependent  variable  is  written  as  a straight-forward  expansion  of  the 
form 

(x , y » t ) = q±(y)  + sqil(x,y,t)  + s2qi2(x,y,t)  + (2.5.1) 


i = 1,2 
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where  s <<  1 is  a dimensionless  small  parameter.  It  may  be  recalled 
here  that  the  steady-state  problem  has  only  one  independent  variable  y. 
Substituting  the  above  expansion  into  the  governing  equations  (2.3.1) 
through  (2.3.8)  it  is  found  that  the  zeroth  order  problem  is  given  by  the 
steady-state  governing  equations  (2.2.1)  through  (2.2.8).  The  first 
order  problem  is  given  by 


Liquid : 


3u  3v 

±±-  + = 0 

3x  3y 

Jn  . - 3un  . i 3pu  , (,2un 

r + °i  — + Vu 

3vn  , - *11  1 3pn  , , fa2vn 

3t  U1  3x  Pl  3y  l[ax2  dy2 

hi , - 9Tn , h rTn , 32th| 

: 1 3x  1 11  P1cpl[3x2  3y2 


(2.5.2) 


,2ul?) 


(2.5.3) 


(2.5.4) 


(2.5.5) 


3u21  9v21 
3T1+  37^=  ° 

- 3u21  1 3p21  . (32“21  , ’Si' 

+ “2  + V21  - - t2  tt-  + "a  ^r-, 


1 3p21 


82v21  , 3Sl 

+ v„  + 


(2.5.6) 


(2.5.7) 


(2.5.8) 


!!2i  + G !!n  + i,v  _ A-fei  + 

3t  U2  3x  2V21  P2c2L2  3y2  J 


(2.5.9) 
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where  primes  denote  derivatives  w.r.t.  y. 

The  next  step  is  to  substitute  the  expansion  (2.3.1)  into  the 
boundary  conditions  (2.4.1)  through  (2.4.14).  It  may  be  recalled 
that  the  interface  boundary  conditions  (2.4.3),  (2.4.4),  (2.4.7), 
C2.4.8),  (2.4.10),  (2.4.13)  and  (2.4.14)  are  applied  at  the  unknown 
interface  y = n(x,t).  Therefore,  consistent  with  the  small  pertur- 
bation approach,  it  is  assumed  that  n(x,t)  is  a small  disturbance 
on  the  steady-state  value  q = 0 and  that  it  can  be  expanded  as 


n(x,t)  = sn1(x,t)  + s1 2n2(x,t)  + — (2.5.10) 

Since  q(x,t)  is  small,  it  is  permissible  to  transfer  the  above- 
mentioned  boundary  conditions  from  the  unknown  interface  y = q(x,t)  to 
the  known  steady-state  value  y = 0,  again  consistent  with  the  small 
perturbation  approach.  This  is  accomplished  through  a Taylor  series 
expansion  about  y = 0, 


q^x.y.t)  = qi(x,0,t)  + y + 


o ' + 


i - 1,2 


=0  'y=o' 


Or  for  y = n 


q1(x,n,t)  = qi(x,o,t)  + jj-  n 


1 + 1 ni  + 

— n + tt  + 


~ 2 2 . 

=0  3y  y=0 


— i = 1,2  (2.5.11) 


where  n(x,t)  is  given  by  (2.5.10). 


r 
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Substituting  Eqs.  (2.5.10)  and  (2.5.11)  into  the  boundary 
conditions  (2.4.1)  through  (2.4.14)  and  collecting  coefficients  of  s, 
it  can  be  verified  that  the  zeroth  order  problem  is  given  by  the 
steady-3tate  boundary  conditions  (2.2.7)  through  (2.2.16).  The  first 
order  problem  is  given  by 

1.  u1;L(x,y,t)  = 0,  y = -h  (2.5.1 


u21(x,y,t)  = 0,  y = 6 


u2i " un  = (ui  ■ u2)ni  * y = 0 


(2.5.12) 

(2.5.13) 

(2.5.14) 


8u2i  ^v21  ^u11  3vh 

4-  W »T  + ^r  + = <“l“l - U2U2  >V  y =0(2.5.15) 


T11(x,y,t)  = 


-(x,y,t)  = 


(2.5.16) 


6. 

T21(x,y,t)  = 0,  y = <$ 

(2.5.17) 

7. 

k 3Tzl 
23y 

- “Uy11  ' <kl T1  - Vi'  >h‘  y ' ° 

(2.5.18) 

8. 

T2i " Tn  = (Ti  “ y = 0 

(2.5.19) 

9. 

P21(x,y,t)  = 0,  y = 6 

(2.5.20) 

10. 

rnta  " (p21 

3v21  3viJ 

' Pll>  + ^2  - Pi)nl  “ 2p  3y  U1  3y  J 

(2.5.21) 

y = 0 
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11. 

vn  (x,y*  t) = °» 

y = -h 

(2.5.22) 

12. 

v21(x,y,t)  = 0, 

y = 6 

(2.5.23) 

13. 

nlt  + nlx“l  - V11 

H 

O 

«• 

II 

o 

(2.5.24) 

14. 

nlt  + nlxU2  " V21 

0,  y = 

(2.5.25) 

An  inspection  of  Eqs.  (2.5.2)  through  (2.5.9)  and  (2.5.12) 
through  (2.5.25)  makes  it  clear  that  the  small  perturbation  assumption 
results  in  linearization  of  the  unsteady  problem. 

2.6  Travelling  Wave  Solution  of  Small  Perturbation  Equations 

The  small  perturbation  equations  (2.5.2)  - (2.5.9)  and  the 
boundary  conditions  (2.5.12)  - (2.5.25)  exhibit  two  important  properties, 
(i)  linearity  and  (ii)  the  coefficients  of  the  unknowns  and  their 
derivatives  are  either  constants  or  functions  of  y at  most.  The  latter 
property  is  a consequence  of  the  steady-state  parallel  flow  assumption. 
These  observations  suggest  a solution  by  separation  of  variables  of 
the  type 

qi;L(x,y,t)  = q^y)  f(x,t),  i = 1,2  (2.6.1) 

A convenient  functional  form  of  f is  chosen  in  what  follows. 

In  the  boundary  conditions  of  Sec.  2.5  there  appear  terms  in 
and  its  derivatives.  Now  q^(x,t)  is  an  unknown  and  its  suitable 
form  mu3t  be  assumed  subject  to  the  condition  of  boundedness  w.r.t. 
x tn.J  t (see  Sec.  2.4).  Suppose  that  initially  the  interface  is 


disturbed  such  that  it  is  sinusoidal  in  form,  thus 


n (x , 0)  = sheikx 


(2.6.2) 


where  s<<1  is  the  dimensionless  small  parameter  encountered  earlier 
in  Sec.  2.5  and  h is  the  liquid  depth.  Eq.  (2.6.2)  implies  that 
ri/h  <<  1 i.e.  the  amplitude  of  the  sinusoidal  disturbance  on  the 
interface  is  much  smaller  compared  to  the  liquid  depth,  k is  the  wave 
number  of  the  disturbance  given  by 

k = — (2.6.3) 


where  X is  the  disturbance  wavelength. 

Now,  Eq.  (2.6.2)  suggests  a travelling  waveform  for  n(x,t)  - 


n(x,t) 


, i(kx-ujt) 
she 


she 


ik(x-ct) 


(2.6.4) 


where 


(2.6.5) 


is  the  speed  of  propagation  of  the  wave  disturbance  and  w is  the 
frequency. 

It  is  clear  from  the  comparison  of  Eqs.  (2.6.4)  and  (2.5.10)  that 

n1(x,t)  = hei(kx^t}  (2.6.6) 

It  is  now  obvious  that  in  Eq.  (2.6.1) 

f (x,c)  - 


resulting  in  the  solution  form 
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q110*.y,t)  = qi(y)ei(kx"“t)  i = 1,2  (2.6.7) 

It  may  be  mentioned  at  this  point  that  the  work,  of  Secs.  2.5  and 
2.6  is  equivalent  to  assuming  the  following  form  of  solution  from  the 
outset 

q^ (x , y , t ) = q^(y)  + sq^yie1^  + 0(s2)  i = 1,2  (2.6.8) 

Substitution  of  Eqs.  (2.6.6)  and  (2.6.7)  into  the  governing 
equations  (2.5.2)  - (2.5.9)  gives  the  result 
Liquid : 


1. 

iku^  + vj  = 0 

(2.6.9) 

2. 

ikpl 

-iwu,  + iku-.u,  + v.ui  - - + v.  (u.  k2u. ) 

1 1111  ii  i 

(2.6.10) 

3. 

pj 

-iwv..  + ikv.u..  = - — + v1  (vV  - k2v. ) 
1 11  l l i 

(2.6.11) 

4. 

k, 

-iwT  + ikT-u,  + v.T'  = (TV  - k2T, ) 

1 11  11  p.c  , 1 1 

1 pi 

(2.6.12) 

Gas : 

5. 

iku2  + v^  = 0 

(2.6.13) 

6. 

ikp2 

-iuju„  + iku„u„  + v.ui  - - + v„(u_  ~ k2u7) 

2 2 2 22  P2  22  2 

(2.6.14) 

7. 

p2 

-iwv„  + ikv„u„  = + v0(v"  - k^v_) 

2 2 2 2 2 2 

(2.6.15) 

8. 

-iwT2  + ikT2u2  + v2T2  - p c (T2  k 7 2^ 

(2.6.16) 

2 p2 
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Similarly,  substitution  of  Eqs.  (2.6.6)  and  (2.6.7)  into  the 
boundary  conditions  (2.5.12)  - (2.5.25)  gives  the  result 


1. 

^(y)  = 0,  y = -h 

(2.6.17) 

2. 

u2(y)  = 0,  y = 6 

(2.6.18) 

3. 

u2  “ ui  = (u[  - uph,  y = 0 

(2.6.19) 

4. 

U2(u2  + ikv2)  “ pl^u|  + ikvp  = (piui  “ U2uph*  y = 

0 (2.6.20) 

5. 

T2(y)  = 0,  y - -h 

or 

(2.6.21) 

T{Cy)  = 0,  y = -h 

6. 

T2(y)  = 0,  y = 6 

(2.6.22) 

7. 

k2T2  " kiTi  = (V'i  - k2*2)h»  y = 0 

(2.6.23) 

8. 

t2  ~ Ti  = (T'x  ~ T'2)h,  y = 0 

(2.6.24) 

9. 

P2(y)  = 0,  y = 6 

(2.6.25) 

10. 

rk2h  + (p2  - p'^)h  + (p2  - px)  - 2(p2v'2  - PjV^), 

(2.6.26) 

y = 0 

11. 

vx(y)  = 0,  y = -h 

(2.6.27) 

12. 

v2(y)  = 0,  y = 6 

(2.6.28) 

13. 

ikh(u1  - ^)  * vx  - 0 

(2.6.29) 
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14. 


ikh(u2  - £ ) - v2  = ° 


(2.6.  30) 


Thus  the  original  unsteady  problem  in  three  independent  variables 
x,  y and  t has  been  reduced  to  a problem  in  one  independent  variable  y. 
The  total  order  of  the  system  of  Eqs.  (2.6.9)  - (2.6.16)  is  14  and 
there  are  14  boundary  conditions  (2.6.17)  - (2.6.30),  thus  the 
resulting  problem  is  mathematically  well-posed. 


2.7  Reduction  to  Four  Dependent  Variables 

There  are  eight  unknowns  u^,  u2>  p^,  p2>  T^,  T2>  v^,  v2  in  the 
governing  equations  and  boundary  conditions  of  the  previous  section. 
Following  the  standard  procedure  in  boundary  layer  stability  theory 
p^  is  eliminated  between  Eqs.  (2.6.10)  and  (2.6.11)  by  differen- 
tiation and  u^  is  eliminated  from  the  resulting  equation  through  Eq . 
(2.6.9).  This  manipulation  results  in  a single  ordinary  differential 
equation  in  v^.  Identical  operations  on  the  gas  side  equations 
(2.6.13)  - (2.6.15)  yield  an  ordinary  differential  equation  in 
Thus  the  system  of  governing  equation  reduces  to  the  following  — 


Liquid : 


1.  v*V-  2k2v1"  + k4v1  » [(v;  - k2v1)(u1  - w/k)  - u^'vj 

k 


-iuT.  + ikT.  u.  + v.T  ’ = — - — (T ’’  - k2T,) 
1 J-l  1 1 p -c  , 1 1 


lpl 

3.  vf  - 2k2v2  + kl|v2  = —■  [(v2  - k2v2)(u2  - w/k)  - u£v2J 


(2.7.1) 

(2.7.2) 

(2.7.3) 
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4. 


-iwT2  + ikT2u2  + v2T2  = 


k? 

— CT"  - k2T  ) 

P2Cp2  2 2 


(2.7.4) 


I 


Eqs.  (2.7.1)  and  (2.7.3)  are  the  well-known  Orr-Sommerf eld  equations. 
Notice  that  the  order  of  the  system  has  been  reduced  from  14  to  12  due 
to  elimination  of  and  p2> 

In  boundary  conditions  (2.6.17)  - (2.6.30)  u^,  u2 , p^  and  p2  are 
eliminated  using  Eqs.  (2.6.9),  (2.6.10),  (2.6.13)  and  (2.6.14).  For 
instance,  solving  for  u-^  and  u2  from  Eqs.  (2.6.9)  and  (2.6.13), 


U1  = " ik 


u2  ' ik 


Solving  (2.6.10)  for  p1  and  using  (2.7.5) 

Pn  . pivi 

7Cv’"  - kV  +~nr(ui  ~ u/k) 


Pi 


P1V1U1 

ik 


(2.7.5) 

(2.7.6) 


(2.7.7) 


Similarly,  from  (2.6.13)  and  (2.6.14) 


-2  - + 


P2V2  - 

4iT(u2  ~ u/k) 


P2V2U2 

ik 


(2.7.8) 


Substituting  Eqs.  (2.7.5)  - (2.7.8)  into  the  boundary 
conditions  (2.6.17)  - (2.6.30)  the  results  are 


1. 

o 

II 

> 

2. 

^2  = °» 

= -h 


= 6 


(2.7.9) 

(2.7.10) 
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3. 

vi  ~ v2  = ikh(ux  - upf  y = 0 

(2.7.11) 

4. 

Px(v^  + k^)  - h2Cv^  + k2v2)  = ikh(y ^u^1  - h2u'p  * 

y = 0(2.7.12) 

5. 

Tj_  = 0,  y = -h 

or 

= 0,  y = -h 

(2.7.13) 

6. 

T2  = 0,  y = 6 

(2.7.14) 

7. 

k2T2  " Vl  = (kl*l  ~ k2^2  y = ° 

(2.7.15) 

8. 

t2  - Ti  = (T|  - T2)h,  y = 0 

(2.7.16) 

9. 

Tk2h  + h(p2  - P^>  + — [m2(v2"  - k2v2)  - Mx(vx' 

k2 

- k2vp] 

+ { [p2^2"2  " v2(u2  ‘ W/k)}  " pl{vlul  " V1(G1  ~ 

w/k) }] 

- 2 (p2v2  - Mjvj)  = 0 at  y = 0 

(2.7.17) 

10. 

v2  = 0,  y = -h 

(2.7.18) 

11. 

v2  = 0,  y = 6 

(2.7.19) 

12. 

v^  - ikh(u^  - u>/k)  =0,  y = 0 

(2.7.20) 

13. 

v 2 - ikh(u2  - w/k)  = 0,  y = 0 

(2.7.21) 

Several  important  observations  can  be  made  at  this  stage, 
(i)  The  number  of  boundary  conditions  has  gone  down  from 
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14  to  13.  This  is  because  p2  no  longer  appears  as  an  unknown  in 
Eqs.  (2.7.1)  - (2.7.4),  consequently  the  boundary  condition  (2.6.25) 
is  superfluous.  Thus  for  a 12th  order  system  there  are  13  boundary 
conditions. 

(ii)  The  governing  equations  (2.7.1)  - (2.7.4)  are  homogeneous 
in  v^,  T^,  and  T£.  The  boundary  conditions  (2.7.9)  - (2.7.21), 
however,  are  not  all  homogeneous;  e.g.  Eqs.  (2.7.11),  (2.7.17), 
(2.7.20)  and  (2.7.21).  This  fact  is  very  significant. 

(iii)  It  is  possible  to  solve  for  v-^  and  v 2 independent  of 
T^  and  T2*  For  instance,  Eqs.  (2.7.1)  and  (2.7.3)  can  be  solved 
subject  to  the  nine  boundary  conditions  (2.7.9)  - (2.7.12)  and  (2.7.17) 

t-  (2.7.21).  This  means  that  the  energy  equation  is  decoupled  from 
the  equations  of  motion.  It  will  be  shown  in  Chapter  III  that  such 
decoupling  is  not  possible  when  there  is  mass  transfer  across  the 
interface.  Once  a solution  for  v^  and  V2  is  obtained  u^,  and  p^,  P2 
can  be  obtained  from  Eqs.  (2.7.5)  - (2.7.8),  if  necessary. 

(iv)  Suppose  Eqs.  (2.7.1)  and  (2.7.3)  are  solved  subject  to 
the  eight  boundary  conditions  (2.7.9)  - (2.7.12)  and  (2.7.17)  - 
(2.7.20)  to  obtain  a general  solution  of  the  form  — 


v1  "*  8^  I y ; P 2 » P 2 9 * ^2  * ^ 9 ^ * ^e  ’ ^ 8^  t P ^ » P 2 * P ^ > ^2 9 ^ ~ 


at  the  interface 


v2  = g2ty;p1,P2,P1,U2,h,6,r,ue,u,k]  = g^p^p^bp^.h,  6,  r,ue,w  ,k] 


at  the  interface 
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r 


All  the  arguments  of  and  g^  except  « and  k are  fluid  properties  and 
remnants  of  upstream  history  (h,6  and  ug)  which  are  known  for  a given 
problem.  However,  w and  k are  not  independent  — they  are  connected 
through  the  last  boundary  condition,  Eq.  (2.7.21).  Thus  this 
equation  is  like  a characteristic  or  frequency  equation  and  in 
dimensional  form  it  is  written  as 

GCP1,P2»M1,U2»h»(S»r»ue»a,,k)  = 0 (2.7.22) 

In  the  present  work,  G will  be  called  the  characteristic  function. 

It  is  clear  from  (2.7.22)  that  given  k,  u is  uniquely  determined 
(there  may  be  more  than  one  value  of  oj)  for  a given  set  of  parameters 
and  vice  versa.  In  this  sense  the  present  problem  is  an  eigenvalue 
problem. 

(v)  Eq.  (2.7.22)  suggests  the  following  method  for  investi- 
gating the  stability  of  the  interface.  Suppose  that  it  is  desired  to 
know  whether  the  interface  is  stable  with  respect  to  a disturbance  of 
wavelength  Given  \ (and  hence  k)  it  is  possible,  in  principle,  to 
determine  a set  of  values  of  <*>.  These  values  of  w will  be,  in  general, 
complex.  Let  “ = + iw^,  Then  the  equation  of  the  interface 

(2.6.4)  becomes 

. . , u)4t  i(kx-<ort)  co  n ton 

n(x,t)  = she  1 e (2.7.23) 


Hence 
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If  w > 0 interface  amplitude  will  grow  exponentially  with  time, 
i.e.,  unstable  interface 

If  < 0 interface  amplitude  will  decay  exponentially  with  time, 
i.e.,  stable  interface 

If  * 0 interface  amplitude  remains  constant  , i.e.,  stable 
interface 

Let  Eq.  (2.7.23)  be  rewritten  as 

n(x,t)  = Shekciteik(x-crt)  (2.7.24) 

where 

c = u/k  (2.7.25) 

Thus  c has  the  dimension  of  speed  and  it  is  referred  to  as  the  phase 
r 

speed,  c^  appears  in  the  amplitude  term  and  is  calTed  the  amplification 
factor. 

2.8  Non-dimensionalization  of  the  Eigenvalue  Problem 

The  vertical  co-ordinates  in  the  liquid  and  gas  are  non-dimension- 
alized  with  respect  to  the  liquid  depth  and  the  boundary  layer 
thickness  respectively.  Thus 

5 = J (2.8.1) 

n = f (2.8.2) 

The  velocities  on  the  liquid  side  are  made  dimensionless  relative  to 

if 


. 

: 


the  interface  velocity  u 


in  Eq.  (2.2.25)  and  the  gas  side  velocities 
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are  made  dimensionless  with  respect  to  the  edge  velocity  ug.  Hence 
the  steady-state  velocity  profiles  now  have  the  form 

Liquid  velocity  profile: 

fljCO  -1  + C • 1 < t 0 (2.8.3) 

Gas  velocity  profile: 

G_(n)  = ^ 0 < n < 1 (2.8.4) 

2 1+E  p — — 


The  interface  velocity,  non-dimensionalized  with  respect  to  boundary 
layer  edge  velocity,  is  given  by 


u = 


(2.8.5) 


where  the  non-dimensional  thickness  and  viscosity  ratios  e and  p 
are  given  by 


(2.8.6) 

(2.8.7) 


The  next  step  is  to  non-dimensionalize  the  governing  equations 
(2.7.1)  and  (2.7.3)  and  the  boundary  conditions  (2.7.9)  - (2.7.12) 
and  (2.7.17)  - (2.7.21).  To  this  end  the  following  quantities  are 
introduced 


(2.8.8) 
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i 2 

*2  " T 

e 


(2. 


Then  using  Eq.  (2.8.1)  and  (2.8.2)  the  Orr-Sommerf eld  equations 
(2.7.1)  and  (2.7.3)  become  (for  a linear  velocity  profile) 


*1V  ~ 2oti^i  + “i^i  = ialRl(lJ;l  "al^l)  (<il  " cl) 


(2. 


and 


i>2  - ^2^2  + a2^2  = ia2R2^2  “ a2^2^  ^U2  ~ c2^ 


(2. 


where  primes  and  dots  denote  differentiation  w.r.t.  £ and  q 
respectively. 

dimensionless  disturbance  wave  number  in  the  liquid 


liquid  Reynolds  nun  her 


cx-^  = kh 


h 


dimensionless  phase  speed  in  the  liquid 


c 


1 


to/k 

Uif 


dimensionless  disturbance  wave  number  in  the  gas 


(2 


(2 


(2 


8.9) 


8.10) 


8.11) 


8.12) 


8.13) 


.8.14) 


a 2 = k6 


(2.8.15) 
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gas  Reynolds  number 


2 v. 


(2.8.16) 


dimensionless  phase  speed  in  the  gas 


c„  = — 
2 u 


(2.8.17) 


The  following  relationships  exist  between  *a2’cl ,C2’  and  R-^.R.^ 


= ea2 


c2  UC1 


r - R 

R1  " - R2 


(2.8.18) 


(2.8.19) 


(2.8.20) 


The  boundary  conditions  (2.7.9)  - (2.7.12)  and  (2.7.17)  - (2.7.21) 


assume  the  following  form  for  a linear  velocity  profile. 


¥[(0  = 0,  t = -1 


2.  ^2(n)  =o,  n = l 

3.  ul^(a  - iOjUj]  = el<l'2(n)  - ia^]  at  £ = 0,  n = 0 

4.  uI^(S)  + a^1(C)]  = ITe2 1^2  (n)  + a^2(n)]  at  £ = n = 0 


(2.8.21) 


(2.8.22) 


(2.8.23) 


(2.8.24) 


».  — — | $2 (n)  - ot2ii«2 (n)  J 


P cfR, 


j ij^U)  - a^x(0 


+ £-  |<wn)  - (G2(n) " c2)^n))'  - wi(5)  - ci>^i<^>} 


J 
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- {e^2Cn)  - u*la)}  ' " |rKw2  + p> 


6. 

7. 

8. 

9. 

where, 

Weber  number 


^(O  = o,  e = -1 


^2^)  = 0,  n = l 


- ia1(fl1(C)  - Cj^)  =0  at  5 = 0 


^2(n)  - ia1(Q2(n)  - c2)  = 0 at  n = o 


“ * 

lpnJ 


(2.8.25) 

(2.8.26) 

(2.8.27) 

(2.8.28) 
(2.8.29) 


p»lf2h 


(2.8.30) 


Froude  number  F = u±f/ 


(2.8.31) 
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3.1  Simplifying  assumptions 

The  following  simplifying  assumptions  were  made  for  the  mass 
transfer  problem  in  order  to  obtain  a mathematically  tractable  model. 

(i)  As  shown  in  Fig.  lb  the  liquid  is  injected  at  y = -h.  At  the 
interface  the  liquid  vaporizes  and  the  vapor  is  entrained  into  the 
gas  boundary  layer  through  convection  and  diffusion.  It  is  assumed 
that  under  steady-state  conditions  the  liquid  injection  rate  exactly 
balances  the  loss  of  liquid  species  at  the  interface.  This  assures 
that  a liquid  layer  of  constant  depth  ' h ' is  maintained. 

(ii)  The  liquid  layer  thickness  h and  boundary  layer  thickness  6 
are  assumed  to  be  prescribed  by  a suitable  upstream  solution. 

(iii)  The  liquid  and  gas  motion  are  assumed  laminar  (or  quasi-laminar) 
and  two  dimensional  both  for  the  steady-state  and  the  unsteady  problem. 
In  the  latter  case  only  two  dimensional  disturbances  are  considered. 

(iv)  The  gas  vapor  mixture  has  constant  properties  such  as  density, 
viscosity,  thermal  conductivity  and  specific  heat.  In  the  case  of 
small  rates  of  mass  transfer,  which  is  the  concern  of  this  work,  these 
properties  have  nearly  the  same  values  as  the  gas  alone. 

(v)  The  Prandtl  and  Lewis  numbers  for  the  gas  mixture  are  unity  in 
both  the  steady-state  and  unsteady  cases.  Thus  the  velocity,  temper- 
ature and  concentration  boundary  layers  have  the  same  thickness  6. 

3 . 2 The  Steady-State  Problem 


fhe  steady-state  or  the  mean  flow  is  assumed  to  be  incompressible 
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and  parallel  with  uniform  injection  rate  at  the  wall  and  uniform 
evaporative  mass  transfer  at  the  interface  (Fig.  lb).  The  mass 
transfer  rate  is  assumed  to  be  consistent  with  the  thermodynamic 
conditions  and  its  determination  is  described  in  Sec.  3.3.  The 
governing  equations  are 
Liquid : 

1.  Continuity 


dy 


0 


(3.2.1) 


2 .  x-momentum 


(3.2.2) 


3.  y-momentum 

dpj 

dy“  = ' Pl8 


(3.2.3) 


4.  Energy 


k1  d2T 


1 


(3.2.4) 


Gas-vapor  boundary  layer: 

5.  Continuity 


0 


(3.2.5) 
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6 . x-momentum 


(3.2.6) 


7 . y-momentum 


dp, 

i 

dy 


= 0 


(3.2.7) 


8.  Energy  (neglecting  viscous  dissipation  and  pressure 
gradient  terms) 

dT„  k„  d2T„ 


2 dy  p„c 


2 p2  dy^ 


n.  2.8) 


9.  Species  continuity 


v,  & - D 

2 d?  dy2 


(3.2.9) 


where  x is  the  mass  fraction  of  vapor  in  the  gas  boundary  layer . 

The  total  order  of  the  system  of  Equations  (3.2.1)  - (3.2.9)  is  14  and 
an  equal  number  of  boundary  conditions  must  be  provided.  These 
conditions  are  listed  below  and  their  explanation  thereafter. 

1.  No  slip  at  the  wall 

u^(— h)  = 0 (3.2.10) 

2.  Boundary  layer  edge  condition  on  the  velocity 


r 
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u2(S)  = ue 


(3.2.11) 


3.  No  slip  in  the  tangential  velocity  at  the  interface 


ux(0)  = u2(0) 


4.  Balance  of  shear  stresses  at  the  interface 


(3.2.12) 


^1 

*1  dy 


du, 

i 

y=0  = V2  dy 


y=0 


(3.2.13) 


5.  Constant  temperature  or  adiabatic  wall 

dT. 


Tx(-h)  = Tw  or  w 


y=-h 


(3.2.14) 


6.  Boundary  layer  edge  condition  on  temperature 


T2<4>  ' Te 


7.  Energy  balance  at  the  interface 

Heat  transferred  from  the  gas  to  the  interface  is  partly 
conducted  through  the  liquid  and  the  remainder  is  spent  in 
vaporizing  the  liquid. 


dT  dT 

ki  dT  ' k2  W + piV  “ y ' 0 


where  l is  the  latent  heat  of  vaporization  of  liquid. 


(3.2.15) 


(3.2.16) 
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8.  No  jump  in  temperature  at  the  interface 


Tx(0)  = T2(0) 


9.  Boundary  layer  edge  condition  on  pressure 


?2(6)  = Pe 


10.  Balance  of  normal  stresses  at  the  interface 


£>i  + = p2  + p2v2  at  y = 0 


(3.2.17) 


(3.2.18) 


(3.2.19) 


11.  Specified  injection  velocity  at  the  wall 


v^-h)  = m/p^ 


(3.2.20) 


where  m is  liquid  mass  transfer  rate  (mass  injected  per  unit  time 
per  unit  area)  at  the  wall. 

12.  Global  mass  balance  at  the  interface 


plvl  = P2V2  at  >'  = 0 


13.  Balance  of  liquid  species  across  the  interface 


P1V1  = P2V2*  " p2°  d$ 


14.  Boundary  layer  edge  condition  on  vapor  concentration 


(3.2.21) 


(3.2.22) 


X(6)  = 0 (3.2.23) 

Note  that  i-  Eqs.  (3.2.9)  and  (3.2.22),  D = v since  = Le2  = 1 





The  first  ten  boundary  conditions  (3.2.10)  - (3.2.19)  are  modi- 
fied forms  of  Eqs.  (2.2.7)  - (2.2.16)  to  account  for  mass  transfer. 


The  value  of  mass  flux  in  at  the  wall  is  determined  by  the  thermo- 
dynamic conditions  of  the  problem  (Sec.  3.3).  Eq.  (3.2.21)  states 
that  the  mass  flux  of  liquid  reaching  the  interface  is  balanced  by 
the  gas-vapor  mass  flux  leaving  the  interface  (note  that  subscript  '2' 
now  denotes  the  gas-vapor  mixture).  Eq.  (3.2.2)  expresses  the  fact 
that  the  mass  flux  of  liquid  species  at  the  interface  is  balanced 
by  the  convection  (p2^2x)  anc*  diffusion  (-p2^9x/3y)  the  vapor 
species.  If  a similar  condition  is  written  for  the  gas  species  at 
the  interface,  assuming  that  the  air  is  insoluble  in  the  liquid  (i.e. 
zero  mass  flux  of  air  at  the  interface), 

0 = p2^2^g  " D2°  d^  (3.2.24) 


where  x Is  the  mass  fraction  of  gas  species.  Since  x = 1 ” X > 
the  above  equation  can  be  written  in  terms  of  x as> 


p2v2  = PpvoX  ~ P,D 


dX 


2 2' 


J2 “ dy 


(3.2.25) 


when  Eq.  (3.2.25)  is  combined  with  (3.2.22)  the  result  is  Eq. 
(3.2.21).  Thus  the  latter  equation  can  be  said  to  express  the 
condition  of  insolubility  of  the  gas  species  in  the  liquid. 
Finally,  the  edge  condition  or  concentration  (3.2.23)  could  well  be 
x(6)  = Xq>  where  Xq  is  some  prescribed  vapor  mass  fraction  in  the 
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inviscid  free  stream.  Thus  x in  the  present  formulation  could  be 
treated  as  a ’reduced'  concentration. 

The  solution  of  the  steady-state  problem  with  mass  transfer 
(as  a function  of  m ) is  — 


v-component  profile  in  liquid 

= m/p^  = constant 

v-component  profile  in  gas-vapor 

v^  = m/p2  = constant 

u-component  profile  in  liquid 

u^  exp(my/p^)  - exp(-mh/p^) 

u exp(m<5/pn)  - exp(-mh/p. ) 

e l -L 


u-component  profile  in  gas-vapor 


u2  exp(my/y2)  - expC-mh/p^) 

u exp(m6/p„)  - exp(-mh/p  ) 

el  J- 


(3.2.26) 


(3.2.27) 


(3.2.28) 


(3.2.29) 


pressure  profile  in  liquid 


p = p + 

K1  re 


1_  1_ 
P2  " o1 


pLgy 


(3.2.30) 


pressure  profile  in  gas-vapor 


P2  Pe  “ 


constant 


(3.2.31) 


H * 


63 


temperature  profile  in  liquid 


■^■|exp(myPr1/M1)  - exp  (-mhPr^p^j  - ^|exp  (rnyP^/^)  - 1 - 


~ (exp(m6Pr2/u2)  ~ 1} 


Y~  ~ c ~Y ~ {exp(myPr1/u1)  - exp (mbP^/p^)  } 


L e pi  e 


T c , 

£ r / • p — i \ i T . P^ 


{exp(m6Pr„/p„>  - 1}  + — {1  - exp (-mhPrj^/u x)  J 

°pl 


constant  temperature  wall 


= 1 - 


— {exp(m6/p„)  - 1}  = const. 

c „T  r L 


(3.2.32) 


p2  e 


adiabatic  wall 


temperature  profile  in  gas 


{exp(myPr2/y2)  - 1}  + -£■ — {1  - exp  (-mbP^/p^)  } + 

‘pl 


{exp(m6Pr2/p2)  - exp (myPr2/y2) } 

”x 

T - r T 1 exp (mbPr1 /u^ ) 

|_  e pi  e J 

{exp(m6Pr2/p2)  - 1}  + (1  - expCmhPr^/u^) } 

Cpl 


constant  temperature  wall 


1 ~ - — Y~  texp(m6/p2)  - exp(my/p2)} 
p2  e 


(3.2.33) 


adiabatic  wall 
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where 


Pr.  = y.c  /k.  = Liquid  Prandtl  Number 
1 1 pi  i 

and 


Pr„  = u0c  0/k„  = Gas  Prandtl  Number  = 1 

2 2 p2  2 


Vapor  mass  fraction  profile 

X = 1 - exp(my/p2)  /exp^/i^) 


(3.2.34) 


The  interface  quantities  are  obtained  from  Eqs.  (3.2.26)  - 
(3.2.34)  by  putting  y = 0.  Thus 


Interface  velocity  — v component 


vif  = ™/pl  y = 0 


(3.2.35) 


m/p  2 y = 0+ 


Interface  velocity  - u component 


uif  1 - exp(-mh/u1) 
ug  exp(m6/p2)  - exp(mh/p^) 


(3.2.36) 


Interface  temperature 


{1  - exp(-mhPr. /y  ) } + 
c ..  11 


if 


ri 

{exp(m6Pr2/y2)  ~ U 

{1  exp ( mhPr1/u1)} 

. e pi  e J 

(exp (mfiP^/l^  ~ 1)  + ll  - exp(-mhPr^/u^) } 


Pi 


(3.2.37) 


W 
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Finally,  one  important  fact  needs  to  be  brought  to  the 
attention  of  the  reader.  In  the  case  of  the  gas  boundary  layer 
the  conditions  u^  * ug,  and  x = 0 were  applied  at  y = 6 

rather  than  at  y + «.  This  was  done  in  order  to  obtain  bounded 
solutions.  Consequently,  these  solutions  have  discontinuous  first 
derivatives  at  the  edge  of  the  boundary  layer. 

3 . 3 Determination  of  Mass  Transfer  Rate  m 

It  was  mentioned  in  the  previous  section  that  the  mass  flux 
m is  determined  by  thermodynamic  conditions.  This  task  is 
accomplished  as  follows. 

It  has  been  assumed  so  far  that  m is  specified  and  that  the 
liquid  depth  h remains  constant.  The  latter  implies  that  whatever 
amount  of  injected  liquid  reaches  the  interface  must  vaporize  and 
then  convect  and  diffuse  into  the  gas.  Now  Eq.  (3.2.34)  shows  that 
at  the  interface  the  vapor  has  a definite  concentration  under  steady- 
state  conditions.  If  the  vapor  alone  were  to  occupy  a unit  volume 
above  the  interface  it  will  be  in  phase  equilibrium  with  the  liquid 
at  the  interface  temperature  and  pressure.  Hence  this  'saturation' 
condition  fixes  the  partial  pressure  of  the  vapor  at  the  interface 
temperature.  The  partial  pressure  of  the  vapor,  in  turn,  determines 
the  interface  concentration.  The  phase  equilibrium  is  expressed  by 
the  Clausius-Clayperon  equation  as 


(3.3.1) 


for  a vapor  point  that  behaves  like  an  ideal  gas.  p^  is  the  partial 
pressure  of  the  vapor  at  temperature  T2*  Thus 


1 


P = Xp- 
*v  t 


(3.3.2) 


by  Dalton's  law  for  an  ideal  gas.  In  (3.3.2)  x is  the  concentration 
and  p is  the  mixture  pressure.  In  (3.3.1)  K is  a constant  and  £ is 
latent  heat  of  vaporization  which  is  assumed  constant,  i is  inde- 
pendently known  from  calorimetric  data.  Combining  the  last  two 
equations 


= Ke  £/R^2 


at  y = 0 


(3.3.3) 


where  is  the  gas-vapor  mixture  temperature  at  the  interface. 
Using  a more  convenient  form  of  (3.3.1)  the  result  is 


R *2  Tref 


at  y = 0 


(3.3.4) 


where  Pref  and  ^ref  are  some  reference  'saturation'  conditions,  i.e., 

T ^ is  the  boiling  point  at  pressure  Pref’^  is  gas  constant  of 

the  vapor. 

Recalling  that  the  steady-state  solutions  were  obtained  in  terms 
of  m,  Eq.  (3.3.4)  becomes 


>(m)p2(m) 


i.  1 


R T2(A>  Trefj 


m can  be  determined,  in  principle,  by  solution  of  (3.3.5).  For 
instance,  in  the  case  of  constant  temperature  wall,  substituting  for 
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p2>  ^2  anc*  X from  Eqs.  (3.2.31),  (3.2.33)  and  (3.2.34)  into  Eq.  (3.3.5) 


RT 


- Un- 


ref 


ref 


£n{l  - exp (~m6/u2)} 


N 

D 


where 


N = exp (mSPr^/v^  ~ 1)  + (1  “ exp  (-mhPr^/y^)  } (3.3.6) 


pl 


D = 


= 

Cpl 


{1  - exp(-mhprl/u1)  } 


+ {exp(m6Pr2/p2) 


- — =-  U - exp  (-mhPr.  /p.  ) } 

C . I 1 -L 

pl  e 


with  Pr2  = 1 

Eq.  (3.3.6)  is  nonlinear  (even  when  simplifications  are  made 
for  small  mass  transfer  rates)  and  is  solved  by  the  Newton-Raphson 
method.  This  procedure  is  described  in  Sec.  5.3. 


3.4  The  Unsteady  Problem 

The  formulation  of  unsteady  problem  follows  closely  the  zero 
mass  transfer  case.  When  the  steady-state  configuration  of  Fig.  1.1b 
is  disturbed  the  resulting  unsteady,  two  dimensional,  incompressible 
motion  with  interface  mass  transfer  is  governed  by  the  equations  below. 


Liquid : 

1.  Continuity 


3u, 


3x 


3v, 


3y 


0 


(3.4.1) 
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8.  Energy  (neglecting  viscous  dissipation  and  pressure  grauient  terms) 


3T2  _ 3T2 

sF  + U2  IT  + V2 


^2  _ 

k2 

32T  32T  ' 

^ -f-  - 

3y 

p2Cp2 

l3x2  3y2 

(3.4.8) 


9.  Continuity  of  vapor  species 


+ u + v 4^  = D 
9t  3x  8y 


fiS  + i5 

3x2  3y2 


(3.4.9) 


This  system  of  equations  is  to  be  solved  subject  to  the 
following  boundary  conditions.  This  development  follows  very 
closely  the  work  in  Sec.  2.4. 

1.  No  slip  at  the  wall 


u^x.y.t)  = 0,  y = -h 

2.  Edge  condition  on  the  u-velocity  component 

u2(x,y,t)  = ue,  y = 6 

3.  No  slip  in  tangential  velocity  at  the  interface  (Fig.  2) 
(same  as  in  zero  mass  transfer  case) 

«2  - “1  " 72}  U at  y = n(x’° 


(3.4.10) 


(3.4.11) 


(3.4.12) 
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4.  Balance  of  shear  stresses  at  the  interface 

(no  modification  over  zero  mass  transfer  case) 


1 - nx 

’3u2 

^2 

2p2nx 

3u2 

3v2' 

i + ni  1,2 

.3y 

■ 3x  . 

1 + n2 

3x 

3y  . 

l - n2 

X 

3u.  3v. 

2Mlnx 

^1 

1 + n2 

X 

3y  3x 

1 + n2 

X 

3x 

9y  . 

5.  Constant  temperature  wall 


T1(3c,y,t)  = Tw,  y = -h 


or  adiabatic  wall 


3T1(x,y,t) 

3y 

6.  Edge  condition  on  temperature 


0,  y = -h 


T2(x,y,t)  = Te,  y = 6 


7.  Energy  balance  at  the  interface 


3T  3T 

on  y ‘ 


(3.4.13) 


(3.4.14) 


(3.4.15) 


where  Vr2  is  gas  velocity  vector  relative  to  the  interface  and 
twai 

the  previous  equation  can  be  written  as 


TT  is  the  outward  unit  normal  to  the  interface.  Since  VR2  = - Vif, 


k2VT2-n  = k1VT1-n  + *p2(V2  - V.f)-n 

If  F(x,y,t)  = 0 represents  the  interface  the  unit  normal  is  given  by 

VF 


n = 


rvFi 


A 


Vh  • ■ 'T1  ' #T  + 7WT  lv2  ' ,F  ■ v< 


A very  delicate  argument  needs  to  be  made  with  regard  to  the  term 

• VF.  At  every  instant  of  time,  the  interface  shape  is  given  by 
a surface  formed  by  those  points  that  have  the  value  F = 0 at  that 
instant.  Thus,  to  an  observer  moving  with  the  interface,  there  is 
no  change  in  value  of  the  function  F.  In  other  words,  the  total  time 
rate  of  change  of  F(x,y,t)  following  a point  on  the  surface  is  zero. 


Hence , 


f + vif  ■ VF  * 0 


The  reader  is  referred  to  Karamcheti  for  a more  complete  discussion 

on  this  point.  Eliminating  the  term  • VF  the  energy  balance 

condition  reads 


k2VT2  • VF  = V^  • VF  + 2p2(V2  • VF  + — ) = 0 


However , 


D„F 

3F  — 3F  — 3F  — 8F  2 

+ V„  * VF  = + u0  i—  + v„  = fTT 

3t  2 3t  2 3x  2 3y  D2 1 


and  hence, 
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D2F 

k2VT2  • VF  = k1VT1  • VF  + lp2  — on  y = n(x,t)  (3.4.16) 

This  equation  reduces  to  Eq.  (3.2.16)  for  steady-state 
conditions  and  to  Eq.  (2.4.7)  in  the  absence  of  mass  transfer. 

8.  No  temperature  jump  at  the  interface 


T1(x,y,t)  = T2(x,y,t)  on  y = n(x,t) 


(3.4.17) 


9.  Edge  condition  on  pressure 


p2(x,y,t)  = p6.  y = 6 


(3.4.18) 


10.  Balance  of  normal  stresses  (momentum  flux)  at  the  interface 
(Fig.  3b) 

This  condition  is  derived  by  applying  Newton's  second  law  to 
the  fluid  crossing  the  interface,  viz. 

Rate  of  change  of  normal  momentum  per  ui  it  area 
= External  stress  in  normal  direction 
Normal  momentum  flux  above  the  interface  = J^(p2^R2  " ^^R2  1 " n 


P2(VR2  ‘ n)2 


Normal  momentum  flux  below  the  interface  = p^(VR^  " 


External  stress  in  normal  direction  = °2  “ °i  _ r 


Hence, 


P2(VR2  ’ *)2  - P1(VR1  ' W)2  = °2  - °1  “ R 
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Following  the  development  of  Eq.  (3.4.16)  it  is  seen  tliat 


V ' n 
R2 


1 V - __  1 V 

JWJ  D2t  and  VR1  ’ n I VF ! Dxt 


Therefore , 


(d  f)2  fn  f)2 

l u2  l _ _ r 

] VFj7  P2  D^t  ~ P1  Dxt  °2  ~ °1  R 


For  F(x,y,t)  = y - n(x,t) 


a + nj)3/2 


^9  ^ /2r  11/2 

2 + f 2 = 1 + n2 

dyj  xj 


Substituting  for  1/R,  | VF | and  2 (from  Appendix  B)  into 

1/2 

Eq.  (A)  and  multiplying  through  by  (1  + n2)  /-  the  following  equation 


is  obtained 


M2 


D-f)2  rn 

1 XX 


2 D„t  MlDlt 


1 1 + n 


- - (p2  - Pl)  (l  + n2) 


. n 2 ^ J 

+ 2n  p„  p.  t — 

x 2 3x  1 3x 


3v2 

+ 2 U2  ~ P1  3jT 

du  3 v.  3u.  SV. 

- 2n  P,  ~ + ~ ~ P,  + T-i 
x 2 3y  3x  1 ( 3 y 3x 


at  y = n(x,t) 


(3.4.19) 
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In  the  absence  of  mass  transfer  Eq.  (3.4.19)  reduces  to  Eq . 
(2.4.10)  and  for  steady-state  conditions  it  reduces  to  (3.2.19). 

11.  Specified  injection  velocity  at  the  wall 

v^(x,y,t)  = m/pj.  y = -h  (3.4.20) 

where  m is  known  from  the  steady-state  solution  of  the  Clausius- 
Clayperon  equation. 

12.  Global  mass  balance  at  the  interface  (Fig.  2) 

This  condition  can  be  expressed  as 


13.  Balance  of  liquid  species  across  the  interface 

This  is  the  condition  that  the  mass  flux  of  liquid  at  the 
interface  is  balanced  by  the  convection  and  diffusion  of  vapor  species 
away  from  the  interface.  Thus 

pl\l  ' " = P2^R2  ' " " d2D  It 

This  equation  reduces  to  Eq.  (3.2.22)  in  the  steady-state  case. 
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Carrying  out  the  usual  substitutions 

1 D1F  - 1 V 1 - 

P1  TW  = P2X  TWT  D2t  P2D  I VF 1"  ' VF 

or 

D F _ D F 

pi  d^t  = P2X  57 ~ p2DVx  • vf  at  y = n(x’° 

combining  with  Eq.  (3.4. 21)  the  final  form  is 

D F 

(1  - x)  + DVX  ‘ VF  = 0 at  y = n(x,t)  (3.4.22) 

Since  it  has  been  assumed  that  Pr^  = Le2  = 1 in  the  unsteady  case  also, 

D = Eq.  (3.4.22). 

14.  Edge  condition  on  vapor  concentration 

x(x,y,t)  =0  at  y = 6 (3.4.23) 

The  boundary  conditions  (3.4.10)  through  (3.4.23)  developed  so 
far  are  based  on  the  same  physical  conditions  as  in  the  steady-state 
case.  The  order  of  the  steady- state  system  of  governing  equations 

(3.2.1)  - (3.2.9)  is  14,  whereas  the  order  of  the  unsteady  system 

(3.4.1)  - (3.4.9)  is  16  with  respect  to  the  variable  y.  As  pointed 

out  in  Sec.  2.4,  the  x and  t dependence  may  be  eliminated  by  assuming 
a travelling  wave  type  of  unsteadiness  and  consequently  one  need  only 
be  concerned  with  boundary  conditions  with  respect  to  y.  The  change 

in  the  order  from  14  to  16  is  due  to  the  appearance  of  second  deriva- 


Cives  w.r.t.  y in  Eqs.  (3.4.3)  and  (3.4.7).  Therefore  two  additional 
boundary  conditions  (on  and  V£)  are  required  for  a well-posed 
formulation.  One  straight-forward  condition,  analogous  to  Eq.  (2.4.12) 
in  the  zero  mass  transfer  case,  is 

15.  ^(x.y.t)  = m/p2 , y = 6 (3.4.24) 

where  m is  the  steady-state  value.  Thus  the  mass  flux  leaving  the 
edge  of  the  boundary  layer  is  assumed  to  remain  unchanged. 

The  last  remaining  condition  is  not  so  obvious.  It  expresses 
the  fact  that  Eq.  (3.4.21)  can  in  fact  be  looked  upon  as  two  boundary 
conditions,  viz. 

D1F  - 

12a.  P1  — = m(x,y,t)  at  y = n(x,t)  (3.4.21a) 

D2F  - 

12b.  P2  ^ = m(x,y,t)  at  y = n(*,t)  (3.4.21b) 

where  m is  the  unsteady  mass  flux  at  the  interface  and  its  determination 
will  be  discussed  later  in  this  section.  Eqs.  (3.4.1)  - (3.4.9)  together 
with  the  boundary  conditions  (3.4.10)  - (3.4.20),  (3.4.21a),  (3.4.21b) 
and  (3.4.22)  - (3.4.24)  form  a well-posed  problem. 

In  the  steady-state  problem  the  mass  transfer  rate  m was  ob- 
tained by  applying  the  condition  of  phase  equilibrium  at  the  interface. 

It  is  now  assumed  that  the  equilibrium  of  phases  prevails  in  the 
unsteady  case  also.  This  requires  that  the  Clausius-Clayperon  equation 
be  satisfied  in  the  unsteady  case.  Writing  Eq.  (3.3.3)  for  the  unsteady 
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problem 

16.  xp2  = Ke-<,/,RT2  at  y = n(x,t)  (3.4.25) 

The  Clausius-Clayperon  equation  thus  determines  the  steady 
(or  mean)  mass  transfer  rate  m for  the  steady-state  problem  and  the 
perturbed  mass  transfer  rate  m in  the  unsteady  case.  With  the  aid 
of  Eq.  (3.4.25)  the  previously  stated  formulation  can  be  modified  as 
follows . 

A well-posed  unsteady  formulation  is  represented  by  governing 
Eqs.  (3.4.1)  - (3.4.9)  and  the  boundary  conditions  (3.4.10)  - (3.4.25). 
It  should  be  noted  that  this  does  not  require  the  use  of  Eqs.  (3.4.21a) 
and  (3.4. 21b)  . 


3 . 5 Travelling  Wave  Solution  of  the  Unsteady  Problem 


The  solution  of  the  unsteady  problem  closely  follows  the  pro- 
cedures in  Secs.  2.5  and  2.6.  The  experience  gained  in  the  solution 
of  the  unsteady  zero  mass  transfer  problem  suggests  the  following 
travelling  wave  solution. 


q (x,y,t)  = q . (y)  + sq  (yje1^*  + 0(s2)  (3.5.1) 

1 1 i=l ,2 


where  s<<l 

The  possibility  of  assuming  the  above  form  of  solution  was 
mentioned  earlier  in  connection  with  Eq.  (2.6.8).  The  interface 
shape  is  assumed  to  be 
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n(x,t)  = she*^*  U)t^  + 0(s2)  (3.5.2) 

The  transfer  of  boundary  conditions  from  the  unknown  interface 
y = n(x,t)  to  the  known  steady-state  position  y = 0 is  accomplished 
by  the  Taylor  series  expansion 

?,  (x,n,t)  = q (x.o.t)  + n + — ^ jr  + — , i=l,2  (3.5.3) 

1 1 3y  y=0  3y2  >=0 

where  n is  given  by  (3.5.2) 

Substituting  Eqs.  (3.5.1)  - (3.5.3)  into  the  governing  Eqs. 
(3.4.1)  - (3.4.9)  and  the  boundary  conditions  (3.4.10)  - (3.4.25)  and 
subtracting  the  zeroth  order  (i.e.  steady-state)  equations,  the  0(s) 
problem  is  (after  considerable  algebra) ; 


Liquid : 

1. 

iku^  + vj  = 0 

(3.5.4) 

2. 

ikPl 

-iu ui  + lkuiui  + viu{  + viui  = - Pl  + vl(ul  - 

k2u:) 

(3.5.5) 

3. 

P{ 

-^i  + ikVi  + Vi  = ~ r-  + vi  (vi  " k2V 

(3.5.6) 

k 

r,  + ikT. u.  + v.T'  +-  v.  T ' = — (TV  - k2T.  ) (3.5.7) 

1 11  11  11  PiC,  1 1 

1 pi 


iku2  + v2  = 0 


(3.5.8) 
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6. 

ikp2 

-lu)u2  + iku2U2  + v2^2  + ^2U2  - ~ + vl^u2  k u2^ 

(3.5.9) 

7. 

p2 

-iuv2  + ikv2u2  + v2v2  = - — + v2^v2  ~ k?v2^ 

(3.5.10) 

8. 

k 

-1wT2  + ikT2u2  + v2T2  + v2X^  - o2Cp2(T2  k T2) 

(3.5.11) 

9. 

- iwy  + ik\u0  + v2x'+  v2x'=  D^x"  “ k^x) 

(3.5.12) 

Boundary  conditions  are 

1. 

u^y)  = 0,  y = -h 

(3.5.13) 

2. 

u2(y)  = 0,  y = 6 

(3.5.14) 

3. 

u2  - ux  = h(u|  - u')  + ikh (v ^ - vp , y = 0 

(3.5.15) 

4. 

P2(u^  + ikv2)  - u1(u{  + ikvp  = (u^"  - P2u2)h’  y = 

0 

(3.5.16) 

5. 

Tx(y)  = 0,  y = 0 

or 

(3.5.17) 

T[(y)  = 0,  y = 0 

6. 

o 

II 

o 

II 

V 

CN 

H 

(3.5.18) 

7. 

k2T2  “ klTl  = £p2[v2  " ikh(“2  “ w'/k'1]  + ^kl^l  ~ k2^")h’ 

y = 0 

(3.5.19) 

8. 

T - T2  - (Tj  - T')h,  y = 0 

(3.5.20) 

9. 


P2 (y)  = o,  y = 6 


(3.5.21) 
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10. 

2m  £ (v2 

-vx)  - 

ikh(u2  - u,)]  = -rk2h  - 

(P2  - PL)  - 

(Pi 

, - 

P[)h 

+ 2 (^2V2  ’ 

- p1v{)>  y = 

0 

(3. 

.5.22) 

11. 

vx(y)  = 0,  y = -h 

(3. 

.5.23) 

12. 

p2[v2  ~ 

ikh(u2 

- u/k)J  = p 1|^v1  - ikh(u1 

- w/k)j  , y 

= 0 

(3. 

.5.24) 

13. 

(1  - 

X)[v2 

- ikh(u2  - w/k)j  - v2x  + 

Dx '=  0,  y = 

0 

(3. 

.5.25) 

14. 

X(y)  = 0,  y = 6 

(3. 

.5.26) 

15. 

v2(y)  = 0,  y = 6 

(3. 

.5.27) 

16. 

Pn  „ W 

^ ~ y = o 

(3, 

.5.28) 

T>2  X RT2 

The  total  order  of  the  system  of  ordinary  differential 
equations  (3.5.4)  - (3.5.12)  is  16  and  there  are  16  boundary 
conditions  (3.5.13)  - (3.5.28). 


3.6  Further  Simplifications 


and  p^  may  be  eliminated  from  Eqs.  (3.5.4)  - (3.5.6)  using 
the  procedure  of  Sec.  2.7  to  obtain  an  equation  in  v^.  Similarly,  an 
equation  in  v2  can  be  derived  by  combining  equations  (3.5.8)  - (3.5.10) 
Now  the  governing  equations  are 


1. 


v?-V-  — v’"  - 2kV;  + “ v . v j + kV 
1 v^  1 1 1 1 1 


= ~ [(V-  - k2v1)(G1  - m/k)  - 


(3.6.1) 
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2. 

-1”1! + “Vi + vi + viTi  • Plcpl(T'i : - kV 

(3.6.2) 

3. 

viv  - -V"-  2k2v"  + V v;  k4v 

2 v2  2 2 v2  2 2 2 

= ^[(V2  ~ k2v2)(u2  ~ u)/k)  ' U2V2^ 

(3.6.3) 

4. 

-i«T,  + ikT  u„  + v.Tl  + v„T ' = (T"  - k2T„) 

2 22  22  22  p„c  „ 2 2 

2 p2 

(3.6.4) 

5. 

-iux  + ikxu2  + v2x'  + v2x'  = D(x"  - k2x) 

(3.6.5) 

A comparison  of  Eq.  (3.6.1)  and  the  original  Orr- Somme rf eld 
equation  (2.7.1)  shews  that  there  are  two  additional  terms  in  the 
former.  These  terms  and  k2v^v|/v^  are  due  to  mass  transfer. 

A similar  observation  can  be  made  by  comparing  Eqs.  (3.6.3)  and  (2.7.2) 


The  variables  u^,  p^,  and  p2  may  be  eliminated  from  the  boundary 


conditions  using 

equations  (3.5.4),  (3.5.5),  (3.5.8) 

and  (3.5.9).  This 

operation  is  the 

same  as  described  in  Sec.  2.7.  The 

resulting  boundary 

conditions  are 

1. 

= 0,  y = -h 

(3.6.6) 

2. 

v^  = 0,  y = 6 

(3.6.7) 

3. 

vj  - v2  = ikh(u|  - u2)  - k2h(v1  - 

v 2 ) , y = 0 (3.6.8) 

4.  Vl(y” 

+ k2v1)  - P2^v2  + k2v2^  = ikh 

- v2u'p , y = 0 

(3.6.9) 

5. 

Tx  = 0,  y = -h 
or 

(3.6.10) 

82 


6. 

t2  = 0, 

y = 6 

(3.6.11) 

7. 

k2T2  " klTl  = ^P2^V2  ~ ikh(^2  ' 

- u/k)J  + (kp^  - k2 

Tph,  y = 0 

(3.6.12) 

8. 

T2  - Ti  = (T*  - 

Tph,  y = 0 

(3.6.13) 

9. 

Tk(i) 2  b + h(p^  - pp  + ^2  [y2(v2"- 

k2vp  - v1  (vp'—  k2vj 

) + m(v'^  - vpl 

+ f [P2*V2“2  _ v2(u2  ” _ P1{V1U1  _ V1<'U1 

- co/k ) >J 

-2(u2V2  - y^vp  + 2m (v 2 - v^)  = 0,  y = C 

) (3.6.14) 

10. 

v1  = 0,  > 

r = -h 

(3.6.15) 

11. 

p2^v2  “ ikh(u2  - w/k)j  = P1^v] 

- ikh(u^  - w/k)J  , ; 

1 = 0 (3.6.16) 

12. 

(1  - X)jv2  - ikh(u2  - w/k)J  - v2x  + Dx ' = < 

3,  y = 0 

(3.6.17) 

13. 

X = o.  > 

f = 6 

(3.6.18) 

14. 

v2  = 0, 

y = 6 

(3.6.19) 

15. 

[y,,t-  k2v'  - JV'l  - ^ 
p2k2L2  2 V22\  ikp2 

i 

r . . l 

V2U2  ' V2(U2  " “/k)J 

tT, 

X 

= , y = 0 

RT2 

(3.6.20) 

The  following  observations  can  now  be  made: 

(i)  The  order  of  the  governing  equation  (3.6.1)  - (3.6.5)  is  14  whereas 


the  number  of  boundary  conditions  is  15. 


The  condition  (3.5.21)  is  no 
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longer  needed  since  P2  has  been  eliminated  as  an  unknown. 

(ii)  The  governing  equations  (3.6.1)  - (3.6.5)  are  homogeneous  in 
Vl’  ^1*  V2*  T2’  X-  Tbe  boundary  conditions  (3.6.6)  - (3.6.20), 
however,  are  not  all  homogeneous.  This  fact  has  important  conse- 
quences as  will  be  shown  in  Chapter  V. 

(iii)  It  appears  that  the  modified  Orr-Sommerfeld  equations  (3.6.1) 
and  (3.6.3)  can  be  solved  independently  of  the  energy  and  concen- 
tration equations.  These  equations,  however,  are  coupled  through 
the  boundary  conditions  and  therefore,  unlike  the  zero  mass  transfer 
problem,  the  present  problem  cannot  be  decoupled. 

(iv)  Eqs . (3.6.1)  - (3.6.5)  can  be  solved  subject  to  14  boundary 
conditions  (3.6.6)  - (3.6.15)  and  (3.6.17)  - (3.6.20)  and  then  Eq. 
(3.6.16)  can  be  used  to  obtain  a characteristic  equation.  This  is 
similar  to  the  method  described  in  Sec.  2.7(iv). 

(v)  The  stability  of  the  interface  is  determined  in  the  same  manner 
as  in  Sec.  2.7(v)  by  solving  for  the  eigenvalue  u>.  Thus 


> 0 


unstable  interface 


o>  < 0 


stable  interface 


u)^  = 0 


neutrally  stable  interface 
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3.7  Non-dimensionalization  of  the  Eigenvalue  Problem 


Non-dimensional  vertical  co-ordinates  £ and  n defined  by 


(3.7.1) 


"-f 


(3.7.2) 


are  introduced.  The  steady-state  profiles  (3.2.26)  - (3.2.34)  are 
made  dimensionless  first,  u^  and  are  made  non-dimensional  w.r.t.  the 
interface  quantities  u^  (Eq.  3.2.36  ) and  (Eq.  3.2.37)  respectively. 
Similarly  u2  and  T2  are  made  dimensionless  w.r.t.  edge  conditions  ug 
and  Tg  respectively. 

Liquid : 


u-velocity  profile 


exp(Rh£)  - exp (-R^) 


Ul(5)  1 - exp(-Rh) 


(3.7.3) 


Temperature  profile  (constant  temperature  wall) 


T 

— w 


c { exp  (R^Pr ^£)  - exp^R^Pr^)}  - c — {exp^Pr-^O  - 1}  + 


‘t 

(expCR^Prj)  - 1}  Ajr-  {exp(RhPri£)  - exp(-RhPri)} 

L e cpl  e J 

r t 

c U - exp (-R^Pr -^) } + {exp(RfiPr2)  - 1>  ^ ^r-{l  “ exP(-RhPrl^ 

p e cpl  e 


(3.7.4) 


w 
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Gas : 


exp(R6n)  - expC-R^) 
U2^n^  exp(Rfi)  - exp(-Rh) 


(3.7.5) 


Temperature  profile  (constant  temperature  wall) 


T2(n)  = 


exptR^Pr^  - 1)  + cp(l  - expt-R^Pr^)}  + 

{exp(R6Pr2)  - exp(R&n)  ) Lp  - ~Y"  U “ exp^R^r^} 
L e pi  e J 

{exp(R(5Pr2)  - 1}  + c"  {l  - exp  (-R^Pr^) } 


(3.7.6) 


Vapor  mass  fraction  profile 


X(n)  = 1 - exp(R^n) /exp(R^) 

(3.7.7) 

where 

Mass  transfer  Reynolds  nunfcer  for  liquid  R^  = mh/p^ 

(3.7.8) 

Mass  transfer  Reynolds  number  for  gas  R^  = mS/p^ 

(3.7.9) 

Liquid  Prandtl  number  Pr^  = pp 

(3.7.10) 

Gas  Prandtl  number  Pr2  = P2Cp2^2  = ^ 

(3.7.11) 

Specific  heat  ratio  c”  = c ~/c  , 

P P2  pi 

(3.7.12) 
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The  interface  velocity  and  temperature,  made  dimensionless  w.r.t. 

boundary  layer  edge  properties,  are 

_ 1 - exp(-Rh> 

U “ exp (Rj ) - exp(-R  ) 


(3.7.13) 


c {1  - cxpC-R^Pr^)}  + (expCR^)  - lip-  _ ^ ~ exP(_RhPrl^ 

— P Le  pi  e J ,3  1 34s 

T " (exp (Rg)  -l}  + c"{l-  expl-R^P^)} 


The  next  step  is  to  non-dimensionalize  the  governing  equations 
(3.6.1)  - (3.6.5)  and  the  boundary  conditions  (3.6.6)  - (3.6.20). 

The  following  dimensionless  terms  are  introduced  for  this  purpose. 


*1  = Vuif 


(3.7.15) 


h = v2/ue 


(3.7.16) 


«1  = Tl/Tif 


(3.7.17) 


82  = VTe 


(3.7.18) 


Then  Eqs.  (3.6.1)  - (3.6.5)  assume  the  form,  for  Pr2  = Le2 


Len  = 1 


1.  ^ - vi""  2aM  + “lVl  + *1*1 

= ic^Rq  - a^1)(u1(5)  - Cq)  - Uq^qj  (3.7.19) 

2.  e"  - Pr^Vl  “ {ai  + ia1Pr1R1(u1(0  - c1)}01  = prqRqVq  (3.7.20) 


3. 


3. 

J2  - R6*2  - 2^2 

+ a2R6^2  + a2^2 

= i«2R2  [^2 

- a^^)  ^u2  “ Cj)  “ u2,i,2 

(3.7.21) 

4. 

02  - ~ (a2  + iot2R2 

(u2  (n)  - c2)}  6 2 = R2T2'*'2 

(3.7.22) 

5. 

'i  - R6X  - (a2  + ia2R2 

(u2(n)  - c2))x  = R2X^2 

(3.7.23) 

where  a^,  a2,  R-^,  R2>  c^>  C2 
(2.8.17)  and  relationships  between 
by  Eqs.  (2.8.18)  - (2.8.20). 

The  boundary  conditions  (3.6 
dimensionless  form 

the  same  definitions  as  in  Eqs.  (2.8.12) 
al’  a2’  cl»  C2  and  R1  R2  are  8iven 

.6)  - (3.6.20)  take  the  following 

1. 

<l>{(5)  = 0 

, e = -i 

(3.7.24) 

2. 

3. 

<i<2(n)  = 0 
u^{(C)  - ia^  + a2  ^J= 

, n = l 

e[hM  - ialG2  + “la2  Rj] 

(3.7.25) 

at  £ = 0,  ri  = 0 

(3.7.26) 

4. 

u |^(0  + a^1(C)j  = yc2^ 

2(n)  + a2iJj2(n)J-  io^GTu^  - 

be  u2) 

at  £ = 0,  n = 0 

(3.7.27) 

5. 

e,  (c)  = 

0,  5 = -1 

(3.7.28) 

6.  0,  (n)  = 0,  n = 1 (3.7.29) 
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7*  ®2(ti)  - el  6ia)  = R2r[*2(n)  ‘ ia1^2CTl)  ^ c2>|+  i[T^  - k‘2*z] 

at  e = 0,  n = o (3.7.30) 

8.  e2(n)  - T01(S)  = TT|  - eT2  at  ( = 0,  n = 0 (3.7.31) 

9.  — i — {^*2 (n)  - a^2(n)}  - = — — {^'"U)  - a^^(C)}  + 

P a?R, 


f2^R7*l(£)  '7i:*2(n)  + i[^*2(n)"2  ' *2(n)(G2(n)  “ 

1 1 / U / 

c2)}- 

= C^i (c)uj[  - ^[(0(^(0  - c1)}j~  -^j^eir<i<2(n)  “ 

+ 

(3.7.32) 

10. 

♦X(C)  = o,  C = -1 

(3.7.33) 

11. 

P^Ol)  “ ia-j^^Cn)  - c2)]=  u[^(C)  - la1(u1(C)  ~ V] 

at  C = 0,  n = 0 

(3.7.34) 

12. 

r2(i  - x(n))^2Cn)  - ia1(u2(n)  - c2)j  - Rgx  + x(n)  = 0 

at  5 = 0,  n = o 

(3.7.35) 

13. 

x(n)  = 0,  n = 1 

(3.7.36) 

14. 

^2(n)  = 0,  n = 1 

(3.7.37) 

15'  ^M*2<n)  " ’ “2*2<n)l‘  T^{“2*2(ri>  - 


/ \ e,(n) 

+ = M— at  c = 0,  n = o 

x(n)  f*(n) 


4*2 (l)  (“2  (h)  * ^2) 


(3.7.38) 


where  in  addition  to  the  non-dimensional  parameters  defined  by  Eqs. 
(2.8.6),  (2.8.7),  (2.8.30),  (2.8.31),  and  (3.7.8)  - (3.7.14)  the 
following  quantities  are  introduced. 


thermal  conductivity  ratio 

k = 

k2/k1 

(3.7.39) 

Euler  number  of  gas 

E = 

Pe/p2Ue 

(3.7.40) 

A = 

*/cplTif 

(3.7.41) 

I = 

a/RTlf 

(3.7.42) 
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4.3.  Mathematical  Statement  of  the  Eigenvalue  Problem 

A concise  mathematical  statement  of  the  problem  described  in 
Sec.  2.8  is  given  below  for  linear  mean  velocity  profiles  described 
by  Eqs.  (2.8.3)  and  (2.8.4). 

Governing  equations : 

1.  t(^V-  2a^  + = ict-jR^i^  - a^Xi^U)  - c^ 

-1  _<  £ < 0 (4.1.1) 

2.  ii>2  - 2 a2  ^2  + = ia2R2^2  “ ^^Hu^n)  - c2) 

0 < n 1 1 (4.1.2) 

where  primes  and  dots  denote  differentiation  w.r.t.  £ and  r) 
respectively. 

Boundary  conditions: 


1. 

<l»j_(-l)  = o 

(4.1.3) 

2. 

^(-l)  = o 

(4.1.4) 

3. 

o 

II 

r—i 

'w' 

CM 

(4.1.5) 

4. 

♦2(1)  = 0 

(4.1.6) 

5. 

u{iJj|(0)  - ic^uj)  = c{^2(0)  - ia^u2 } 

(4.1.7) 

6. 

u{^(0)  + a^1(0)}  = ye2  { 1^(0)  + a^2(0)} 

(4.1.8) 

J 


7.  — - {'^o  CO)  - {*'*(0)  - a2*'(0)} 

l *■  P Jr*  *-  -*-  -*- 


a2R. 


2 2 


^ q2r 

alRl 


TT2  i 


+ a {V2(0)  ‘ (U2(0)  " c2)^2(0)  } pa  {ui*l(0)  * (ul^0) 


{ey^„(0)  - uiJj' (0)  } = - (a2W2  + ) 

A l p l va 


8. 


ip1(0)  - ia1(u1(0)  - 


(4.1.9) 

(4.1.10) 


9. 


*2(0)  - ict1(G2(0)  - c2)  = 0 


(4.1.11) 


where 


Gx(c)  = l + s 

-i  < ^ < 0 

(4.1.12) 

- , \ ey+n 

2(n)  ep+1 

o £ n <_  l 

(4.1.13) 

al 

= ea2 

(4.1.14) 

C2 

= UCj 

(4.1.15) 

R1 

= ueyR^p 

(4.1.16) 

u 

ey 

1+ep 

(4.1.17) 

The  two  Orr-Sommerf eld  equations  (4.1.1)  and  (4.1.2)  are  homo- 
geneous in  and  The  boundary  conditions  (4.1.3)  - (4.1.11), 

however,  are  not  all  homogeneous.  One  of  these  conditions,  say 
(4.1.11),  can  be  used  to  make  Eqs.  (4.1.9)  and  (4.1.10)  homogeneous. 
The  resulting  system  will  then  consist  of  two  homogeneous  fourth 
order  equations  with  eight  boundary  conditions,  i.e.  a legitimate 


i! 


93 


eigenvalue  problem.  In  the  present  work,  however,  the  approach  of 
37 

Bordner  et.al  is  followed.  The  first  eight  boundary  conditions 
(4.1.3)  - (4.1.10)  are  used  to  determine  eight  constants  of  inte- 
gration and  then  (4.1.11)  is  used  to  obtain  the  characteristic 
equation.  This  treatment  was  discussed  earlier  in  Sec.  2.7. 


4.2  Solution  for  a Long  Wavelength  Disturbance 


Consider  a disturbance  on  the  interface  whose  wavelength  is 
much  larger  than  the  liquid  depth  and  the  boundary  layer  thickness. 
Thus 


<<  1 and  (*2  <<  1 


with  a’2./a2  = h/<$  = £ • In  most  problems  of  interest  c itself  is  very 
small  so  must  be  extiemely  small.  Let  and  ij^  be  represented  by 


the  following  straightforward  expansions 


*1  = *10  + al*ll  + “l*12  + — 


*2  = *20  + a2*21  + “2*22  + “ 


(4.2.1) 

(4.2.2) 


Substituting  these  expansions  into  the  governing  equations  and 
boundary  conditions  (4.1.1)  - (4.1.10),  making  use  of  (4.1.14)  and 
equating  coefficients  of  equal  powers  of  a , the  results  are 
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Zeroth  order  problem: 
Governing  equations 


0(a°) 


,1V  n 
*10  = 0 


*20  = 0 


Boundary  conditions 


1. 

*10<-1>  - 0 

2. 

*;0<-i>  - o 

3. 

*20(1)  - 0 

4. 

*20(1)  ' 0 

5. 

U^[0(0)  - ei20(0)  = 0 

6. 

u*iq(0)  - ”^e2*20^^  = 0 

7. 

f><°>  - f k 

8. 

*10(0)  - 0 

(4.2.3) 

(4.2.4) 


(4.2.5) 


(4.2.6) 


First  order  problem:  O(a^) 

Governing  equations 


*11  = iRi(lli  ' Cl}  *10 


*21  = iR2<‘U2  ” C2^*20 


(4.2.7) 

(4.2.8) 
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Boundary  conditions 


*u(-i)  = o 


^(-1)  = o 


*21d>  " ° 


*21(1)  = ° 


R^21(0> 


uiJij^O)  - i21(0)  = i(uG|  - eu2) 

^li^0>  " = 0 

-t  <‘i- 

- ie{u2^20(0)  - (u2  - c2) ^2o (0) ^ 
'i'llC0)  = 1(^(0)  - c1) 


The  characteristic  equation  for  is  obtained  from  Eq.  (4.1.11) 
with  the  aid  of  Eq.  (4.1.15) 

u9(0)  . 

C = — + == Ip  (0)  (4.2.K 

1 U ua^  2 

Substituting  the  expansion  (4.2.2)  into  the  above  equation  and  using 


(4.1.14) 


c,  ■ — 
1 u 


u2(0)  ± 


+ %*20(0)  + W0)  + °‘“1> 


w 


The  zeroth  order  problem  is  completely  homogeneous  and  therefore 


has  only  a trivial  general  solution.  Thus 

*10  2 0 (4.2.12) 

*20  ^ 0 (4.2.13) 

The  first  order  governing  equations  (4.2.7)  and  (4.2.8)  have  the 
general  solutions 

*H  = An  + A12£  + A13£2  + AUC3  (4.2.14) 

*2i  = A2i  + A22n  + A23n2  + A24n3  (4.2.15) 

Introducing  Eqs.  (4.2.14)  and  (4.2.15)  into  the  boundary  conditions 
(4.2.9),  using  (4.1.16)  and  solving  for  the  constants  of  integration, 
the  following  results  are  obtained  after  lengthy  algebraic  manipulations, 

A = i(u  (0)  - c ) 
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i{2(l  + -i)  (uu^  - eu2)  - u(u1(0)  - c^  (3  + p + Hp) } 


(4.2.16) 


i|6u(u^(0)  - Cl)(i  - 1)  - (4  + p)(uu|  - eG2)} 


A23  ' 


l|2(uu|  - eu2)  - 311(1^(0)  - c±)  (1  - ^=)} 


ijCuu^  - eu2)  - 2u(u1(0)  - c.^  (1  + -p=)} 


D = 1 + ^=(4  + f) 

ep  e 


To  complete  the  task  of  obtaining  the  eigenvalue  Eqs.  (4.2.13) 
and  (4.2.15)  are  combined  with  Eq.  (4.2.11)  to  get 


u2(°)  i 

C1  = + ^ A21  + ° (al} 


(4.2.17) 


Finally,  substituting  for  u2(0),  u and  from  Eqs.  (4.1.13),  (4.1.17) 
and  (4.2.16)  into  Eq.  (4.2.17),  the  expression  for  the  eigenvalue  is 

£ *U2  + (2ey  + 1)  {-  + 2 (e  + 1)  } 

c. f : — Z 7 T + 0(o  ) (4.2.18) 

ep  |6  + 4e (1  + --y)  + e2p(l  + —ip?) 

The  expression  for  c^  shows  that  this  particular  eigenvalue 
(or  mode)  depends  only  on  the  thickness  ratio  e and  the  viscosity 
ratio  IT ; it  is  independent  of  the  film  Reynolds,  Froude  and  Weber 
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i 

I 


numbers.  Tne  most  important  observation  is  that  is  purely  real 
and  therefore  represents  a neutrally  stable  mode.  Eq.  (4.2.18)  can 
be  simplified  greatly  for  the  case  of  a thin  liquid  layer  with 
e ,"p  < 1 such  that  eu  <<  1 and  e2  <<  1.  This  operation  results  in 


c1  = 1 + 2e  + 0(a1)  (4.2.19) 

This  form  is  more  illuminating  in  that  it  clearly  shows  that  the 
critical  point  always  lies  inside  the  gas  (i.e.  c^  > 1)  for  this 
mode.  This  statement  will  perhaps  become  more  meaningful  in  Sec.  4.6. 

4.3  General  Solution  for  Arbitrary  Disturbance  Wave  Numbers 


Eq.  (4.1.1)  can  be  written  in  the  form 

- a24'1)"-  <x24"  - a2^)  = iajRj.C*"  - a24-1)(G1(C)  - Cj) 

(4.3.1) 

Let 

^-a2^=w1(C)  (4.3.2) 

Then  Eq.  (4.3.1)  becomes 

W1  ” {al  + iaqRi^1^)  - c^Jw.^  = 0 (4.3.3) 

28 

Following  Feldman  , defining  the  transformation, 

a2  + ia1R1(u1(C)  - c^ 

(o1R1G‘f/3 


CX(C)  = - 


(4.3.4) 
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where  is  a constant  for  a linear  velocity  profile,  Eq.  (4.3.3)  reads 


^lhl  = 


where  h^^)  = w^S^)} 


(4.3.5) 

(4.3.6) 


Eq.  (4.3.5)  is  the  well  known  Airy  differential  equation  and  it  is 
regular  everywhere  in  the  complex  4^  - plane.  In  fact,  this  equation 
possesses  an  irregular  singularity  at  infinity. 

Eq.  (4.3.5)  has  the  following  pairs  of  independent  solutions, 


Ai(^1),  Bi(;l) 

Ai(cl)  , Ai(cie27rl/3) 

AiUp,  Ai(?1e_2lTi/3) 

where  Ai  and  Bi  are  called  the  Airy  functions  of  the  first  and  second 
kind  respectively.  Ai  has  the  property  that  it  is  real  for  real 
and  Bi  is  constructed  from  Ai  in  such  a way  that  Bi  is  also  real  for 
real  The  relationships  amongst  the  above  solution  pairs  are  (Ref.  44) 


Bi 


(Z)  = eni/6  Ai(Ze27rl/3)  + e wi/6  Ai (Ze'2ni/3)  (4.3.7) 


Ai  (Z)  + e2ni/3  Ai(Ze27,l/3)  + e_2l,l/3  Ai(Ze  2it1/3)  (4.3.8) 


Ai(Ze±2lTi/3)  =±e±,Tl/3 


(Ai(Z)  ? iBi(Z)} 


(4.3.9) 
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In  the  present  analysis  the  pair  of  solutions 


Ai(Ze2lTi/3),  Ai(Ze_2lTi/3) 


(4.3.10) 


was  chosen  for  convenience  in  the  numerical  integration  encountered 
later.  Thus  the  solution  to  Eq.  (4.3.5)  is 


h1(C1)  = C3Ai(;1E  ) + C4Ai(SlE_) 


(4.3.11) 


where 


E1  = e±2*i/3 


(4.3.12) 


Hence 


From  Eq.  (4.3.2) 


^(O  = C3Ai{41(OE+}  + C4Ai{c1(C)E-}  (4.3.13) 


- oJ*j_  = C3Ai{^1(C)E+}  + C4Ai{Cl(OE-}  (4.3.14) 


The  homogeneous  solution  of  Eq.  (4.3.14)  is 


^1H  = CiexP(ai^)  + C2exp(-a  £) 


(4.3.15) 


and  the  particular  solution,  obtained  by  the  method  of  variation  of 
parameters  (Ref.  42)  is 


C 

■—  I sinh{a3(£  - t)  }Ai (t) E+}dt 
C f 5 

F I sinh{a1(S  - t) }Ai(41(t)E-}dt 


(4.3.16) 
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and  q*  is  such  that 

c2(n*)  - 0 (4.3.21) 

4.4  Application  of  Boundary  Conditions 

In  the  boundary  conditions  (4.1.3)  - (4.1.12)  derivatives  of 
ip^  and  1^2  w.r.t.  5 and  q,  up  to  third  order,  are  required.  Since 
i and  q occur  in  the  limits  of  integration  of  Eqs.  (4.3.17)  and 
(4.3.19)  it  is  necessary  to  use  Leibniz's  rule  of  differentiation 
under  the  integral  sign.  The  results  are  summarized  in  Appendix  C. 

Substituting  Eqs.  (4.3.17),  (4.3.19)  and  the  derivatives  in 
Appendix  C into  the  boundary  conditions  (4.1.3)  - (4.1.11)  and 
performing  the  required  algebraic  simplifications,  the  following 


equations  are  obtained: 

a^exp(-a^)C^  + a1exp(a^)C2  + XqC3  + I2C4  = 0 (4.4.1) 

a1exp(-a1)C1  - a1exp(a1)C2  + = 0 (4.4.2) 

u2exp(a2)C5  + a2exp(-a2)C6  + + I6Cg  = 0 (4.4.3) 

a2exp(a2)C5  - a2exp(-a2)C6  + + IgCg  = 0 (4.4.4) 


UctlCl  “ ualC2  + TrillC3  + “112C4  “ eoi2C5  + ea2C6  " eI15C7  “ eI16C8 

= ia^ (uu|  - eu2)  (4.4.5) 


Iiuta 
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2ua2C1  + 2Ua^C2  + ^■GA±{c1(0)E+}  - 21*0^1^  C3 
+ jliAlUj.WE"}-  Zuojl^  ] C4  - 2ye2a^C5  - 2ye2a2C6 

^-TIe2Ai{?2(0)E+}  + 2Ue2a2I13J  C?  + ^-ye2Al{ C2  (0)E~ } + 2ye2a2I14J  Cg  = 0 


(4.4.6) 


”2  / \ 4UUi 

^r(Gi-  {Gi(0)  - ci}ai)  + i^ 

. -2  / \ 2ua 

~ ^ f“(Ui  + {ui(0)  " cl}alj  " eTTR2 

-?-lTA1'Ki(0)E+kiE+  + ^^I9 

(f  rlGi(°)  - <=!> + %)  hi] 

+ [*^^4l,Ul(0,ritir+^^I“ 

+(tirlui<0>  * 'i1  +ils“)Ii2] 

2)-?js 

[qt  “■  «,»»%■•  - J ■„  ■(  * y 'u]  ‘i 


+ |^(G2  " {G2(0)  " C2  }a2) 


+ lt^(G2  + {G2(0)  " C2}a2j 

, X 


= - t-  fa2W2  + ^ 


(4.4.7) 
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alCl  + alC2  ~ I9C3  “ I10C4  = ~ c1^  (4.4.8) 

a2C5  + a2C6  - h3C7  - I14C8  “ iola2{a2(0)  " c2 } (4-4'9> 

where  the  integrals  1^  through  are  defined  in  Appendix  E. 

Eqs.  (4.4.1)  - (4.4.9)  form  a system  of  eight  linear  algebraic 
equations  of  the  type 

[a(Ci)](C}  - {V(cj)}  (4.4.10) 

where  £a(c^)J  is  the  coefficient  matrix  of  the  left  hand  sides  and 
(V(c^)}  is  the  column  vector  of  the  right  hand  sides.  The  remaining 
equation  (4.4.9)  can  be  written  in  the  form, 


G [cl;*C^cl^]  a2C5  + a2C6  “ I13C7  ~ I14C8  ~ iai°2^u2^  c2^  “ 0 

(4.4.11) 

or  more  compactly 

G [Cl;{C(Cl)}]  = {C(c1)}T{U(c1)>  - i°‘1a2  (u2  (0)  - c2)  = 0 (4.4.12) 


where 


and 


(C(c1)}  = £a(c^)J  1{V(c1). 


(4.4.13) 


{U(c.)}T  - {0  0 0 0 a,  a,  -I,,  -I,.}  (4.4.14) 

1 2 2 13  14 

In  Eqs.  (4.4.11)  and  (4.4.12)  C2  is  given  by  Eq.  (4.1.15),  G is  the 
characteristic  function  of  Sec.  2.7  and  c^  is  the  eigenvalue. 


The 


problem  of  determining  stability  of  the  interface  is  thus  reduced 
to  finding  the  points  in  the  complex  c^  plane  at  which  G vanishes. 

Eq.  (4.4.12)  can  be  expressed  in  the  functional  form 

G (a^,e  ,W^F;c^)  = 0 (4.4.16) 

For  the  present  physical  problem  G is  an  analytic  function  of  the 
above  parameters  and  the  mathematical  problem  reduces  to  finding  the 
zeros  of  an  analytic  function. 

4.5  Outline  of  the  Eigenvalue  Iteration  Procedure 

The  zeros  of  the  characteristic  function  G must  be  determined 
numerically  and  the  major  steps  in  this  procedure  are  listed  belcw. 

(i)  Guess  c^.  A method  for  obtaining  a good  guess  value  is  described 
in  the  next  section. 

(ii)  For  given  values  of  a^,  e,  p,  p , R^,  W,  F and  the  guess  value 
c^,  evaluate  the  integrals  1^  - I ^ using  a suitable  numerical 
procedure.  This  is  described  in  Appendix  G. 

(iii) Generate  the  coefficient  matrix  [A]  and  obtain  its  inverse.  The 
inverse  was  obtained  using  the  routine  MINV  in  the  IBM  Scientific 
Subroutine  Package  after  making  minor  changes  to  handle  complex 
numbers.  The  constants  of  integration  were  conveniently  scaled  to 
avoid  overflows  in  MINV.  In  order  to  check  the  accuracy  of  matrix 
inversion  the  eigenvalue  was  obtained  using  MINV  alone  and  by 
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applying  a first  order  correction  to  [A]  The  difference  between 
the  eigenvalues  obtained  by  these  two  methods  was  negligibly  small. 
Also  compute  the  right  hand  side  column  vector  {V}. 

(iv)  Determine  the  constants  of  integration  i.e.  the  vector  {C}  by 
carrying  out  the  matrix  multiplication  [A]  ^{V}. 

(v)  Compute  G in  Eq.  (4.4.11),  it  should  be  close  to  zero  if  this 
equation  is  satisfied. 

(vi)  If  Eq.  (4.4.11)  is  not  satisfied  calculate  an  improved  value 
of  c^  by  using  a suitable  technique.  In  the  present  work,  the 
Newton-Raphson  method  was  used  for  this  purpose.  This  iterative 
method  requires  computation  of  the  first  derivative  of  G.  w.r.t.  c^. 
Details  of  the  Newton-Raphson  method  are  given  in  Appendix  I. 

(vii)  Compare  successive  values  of  c^  for  convergence  within  a 
prescribed  tolerance  on  real  and  imaginary  parts.  Repeat  steps 
(ii)  - (vi)  until  desired  convergence  is  reached. 

4.6  Gen  ration  of  an  Initial  Guess  for  c^ 


Since  the  Newton-Raphson  method  is  sensitive  to  the  initial  guess 
it  is  necessary  to  know  an  approximate  value  of  c^.  This  value  could  be 
determined  either  on  a mathematical  basis  or  a physical  basis.  The 
small  perturbation  method  of  Sec.  4.1  can  be  used  as  a mathematical 
basis.  However,  it  leads  to  only  one  zero  of  the  characteristic 
function  G.  It  was  found  in  the  present  investigation  that  a number  of 


t 
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zeros  of  G are  possible  (see  Chapter  6 for  details.) . Hence  a more 
reliable  approach  to  obtain  the  guess  is  necessary.  As  far  as  the 
physical  basis  is  concerned  the  works  of  Lock^  and  Landahl'1'^  suggest 
that  the  real  part  of  c^  or  the  phase  speed  should  be  close  to  the  speed 
of  propagation  of  free  surface  waves.  Guessing  only  the  real  part 
accurately,  however,  is  not  sufficient  because  the  imaginary  part  is 
of  significance  as  well.  It  would  be  adequate  to  guess  only  the  real 
part  if  the  investigation  were  to  be  confined  to  neutral  stability 
analysis.  In  the  latter  case  the  imaginary  part  of  c^  is  necessarily 
zero. 


The  above  discussion  shows  that  a simple  clear-cut  solution  is 
apparently  not  possible.  An  attempt  was  made  to  determine  at  least 
the  number  of  zeros  of  G inside  a closed  contour  in  the  complex  c^ 
plane.  This  was  done  using  the  'argument  principle'  which  states  that 
for  a regular  analytic  function  G(Z)  within  a closed  contour  C the 
number  of  zeu-  of  G within  C is  given  by 


r 


(Z) 

(Z) 


dZ 


(A. 6.1) 


where  it  is  assumed  that  G does  not  vanish  on  C. 

Numerical  evaluation  of  the  above  integral  as  the  limit  of  a sum 


I 


i 


proved  to  be  a very  difficult  task.  This  was  due  to  the  fact  that  the 
integrand  in  (4.6.1)  is  highly  oscillatory  and  undergoes  large  changes 
of  magnitude.  Consequently,  either  graphical  or  purely  numerical 
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determination  of  the  number  of  zeros  in  this  manner  would  be  extremely 
difficult  unless  a very  fine  step  size  along  the  contour  is  employed. 
The  latter  choice,  of  course,  leads  to  large  computation  times.  There- 
v fore,  this  course  also  had  to  be  abandoned. 

The  only  choice  available  is  to  determine  the  zeros  of  Re(G)  and 
Im(G)  separately  and  then  to  obtain  their  common  zeros.  This  can  be 
done  graphically  as  follows.  First  a suitable  interval  on  c^  is 
chosen.  Recalling  that  the  phase  speed  rn/k  was  made  dimensionless 
w.r.t.  the  interface  velocity  {Eq.  (2 . 8 . 14)  } , it  is  seen  that  c = 1 
corresponds  to  a critical  point  at  the  interface.  Thus  disturbances 
which  propagate  with  phase  speed  greater  than  the  interface  velocity 
lie  in  the  interval  1_<  c^r  <_  1/u.  Now  Eq.  (4.1.17)  shows  that  for 
thin  liquid  layers  and  typically  small  viscosity  ratios  ep  <<  1 and 
thus  1/u  - e"p  is  large.  Hence  the  gas  side  interval  on  c^  is  very 
large  compared  to  the  liquid  side  (0  <_  c^  1) . The  critical  points 

of  interest,  however,  are  those  that  lie  near  the  interface.  In  other 
words,  relatively  slow  moving  disturbances  near  the  interface  are  of 
interest.  The  reason  being  that  in  a real  physical  situation  slow 
moving  disturbances  are  more  likely  to  be  triggered.  In  conclusion, 
it  is  sufficient  to  consider  a typical  interval  0 c^^  <_  4. 

The  next  step  is  to  choose  a suitable  interval  on  c^ . It  is  not 
possible  in  this  case  to  offer  an  argument  similar  to  the  previous  one. 
Therefore,  0_<|c^.J_^  1 is  chosen  as  a start.  G is  then  calculated  at 
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a number  of  points  inside  the  unit  rectangle  (e.g. , at  intervals  of 

0.1  along  c^r  and  c^)  . Then  Re(G)  is  plotted  against  Re(c^)  with 

Im(c^)  as  a parameter  and  the  points  of  intersection  on  the  x-axis 

are  determined.  Similar  plots  are  made  for  Im(G)  and  its  zeros  are 

also  determined.  It  is  usually  observed  that  as  one  proceeds  from 

c^^  = 0 to  | c^  | = 1 the  trend  of  the  curves  shews  that  beyond  a 

certain  c,  . there  are  no  intersections  on  the  x-axis.  This  fact 
li 

determines  the  upper  limit  on  c^.  Finally,  common  zeros  of  Re(G) 
and  lm(G)  can  be  roughly  determined.  These  approximate  values  of 
c^  serve  as  initial  guesses  for  Newton-Raphson  iteration. 

An  illustrative  example  of  the  above  procedure  is  included  in 
Chapter  VI. 


CHAPTER  V 


SOLUTION  OF  THE  MASS  TRANSFER  PROBLEM 


Ill 


5.1  Linear  Approximation  of  Exponential  Steady-State  Profiles 


The  present  investigation  concerns  itself  with  small  values  of 
mass  transfer  rates.  Hence  it  is  assumed  that 

Rh  « 1 

and  (5.1.1) 

Rr  « 1 

o 

Since  |c|  and  | n | are  always  less  than  unity  the  exponentials  in 
Eqs.  (3.7.3)  and  (3.7.5)  can  be  expanded  in  a Taylor  series  to  give 
(up  to  first  order) 

UjU)  = 1 + £ , -1  < 5 <0  (5.1.2) 

and 

R,  + R,n 

G2(n)  ~ ?'  + R~  ’ 0 - n - 1 (5.1.3) 

n 6 

It  may  be  verified  that  the  linear  profiles  in  the  last  two  equations 
reduce  to  Eqs.  (2.8.3)  and  (2. 8. A)  with  the  help  of  Eqs.  (3.7.8)  and 
(3.7.9).  The  mass  fraction  profile  of  Eq.  (3.7.7)  reduces  to 


x(n)  = R6(l  - n),  o < n < 1 (5.1.4) 

In  order  to  linearize  (3.7.4)  it  needs  to  be  assumed  in  addition  to 
(5.1.1)  that 

R^)Pr1  <<  1 (5.1.5) 

Typically,  Pr^  - 0(10)  and  this  requires  R^  to  be  smaller  than  0.01. 


2 


:<6  + CpVrl 


o < n < 1 


(5.1.7) 


Finally,  linearized  expressions  for  the  interface  quantities  in 
Eqs.  (3.7.13),  (3.7.14)  and  (3.7.7)  are 


R6  + ®h 


Rr  — — + It  Pr,  c 

6 T Ti  1 p c 

R6  + ?PVrl 


kR*] 


x = Rx 


3.2  Mathematical  Statement  of  the  Eigenvalue  Problem 


The  mathematical  statement  of  the  mass  transfer  problem  of  Sec. 
3.7  is  given  below  for  the  case  of  the  linear  steady-state  profiles 


described  in  Sec.  5.1. 
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Governing  equations: 


1. 

*iv_  VT  ■ H*! + aivl + ai^i 

= ia1R1(^  - a2ip1)(u1(C)  - c.^ 

(5.2.1) 

2.  - 

Pr^Sj  - {a2  + ia1Pr1R1(ii1(5) 

- 'i>6i  • "iVl'i 

(5.2.2) 

3. 

*2  " V2  " 2a2^2  + a2R6^2  + a2*2 

= i a2 R2  ( q32  - a2 

;ii»2)(u2(n)  - c2) 

(5.2.3) 

4.  0*2  - 

R(S®2  ~ ^a2  + ia2R2  ^u2  “ c2 

;)}02  = r2t24< 

(5.2.4) 

X - 

R^X  * (a2  + ia2R2(u2(n)  - c2>  )x  = R2^2 

(5.2.5) 

Boundary 

conditions: 

1. 

^(-1)  = 0 

(5.2.6) 

2. 

<l^(-l)  = 0 

(5.2.7) 

3. 

K»2(l)  = 0 

(5.2.8) 

4. 

*2(i)  = 0 

(5.2.9) 

r 1 

R$" 

5 . "u 

♦i(°)  “ iaiGi  + °i  r:J  = £ 

*2(°)  ' ia2G2  + ala2  rT 

(5.2.10) 
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6. 


u[»l<0)  + - TTe2 r *2(0)  + o|t2(0)j  (5.2.11) 

i-[*2(0)  - «|*2(0)1  -£-i-L"(0>  - »2„;<0)1 

L j aT  -J 


2 ^ 

7-  +rtf  *1(0)  - *2(0)  + 1 ^2(0)li2  - (G2(0)  - c2)^2(0)} 


u2  1 


■ {*l(0)ui  - (ul(0)  - 


Wr{^2(°)-^(0)} 


8. 


9. 


*2<0)  - ~ 


r^*i<0)]  ■ ■ f~  <ai"2  + fj  (5-2-12> 

P’^2<>0^  “ ic*i  (u2  (0)  ~ c2)]  =u^1(0)  - ia^u^O)  - c^l  (5 . 2 . 13) 

R2U  - x(0)}[  *2(0)  - ia1(G2(0)  - c2)]  - R^xCO)  + x(0)  = 0 (5.2.14) 


10. 

11. 

12. 

13. 


e1(0)  = o 


e2d)  = o 


e2(0)  - 10,(0)  = TT ^ - ef2 


02(O) 


~ R,TA 

7xei(0)  = 4" 

P L 


^2 (0)  - ia1(u2(0) 


c2} 


(5.2.15) 

(5.2.16) 

(5.2.17) 

(5.2.18) 


14. 


15. 


a2R 
2 2 


-t*2(0)  - 


V2(0) 


X(D  = 0 


a^2(0)}  - hr 


v2^2 


(5.2.19) 
) 


i>2  (°)  ^u2  (0)  - c2) 


x(0) 


0O  (0) 


t2(0) 


(5.2.20) 
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where  0^(0,  u2(n),  x(n)  , ^(n)  and  T2(n)  are  given  by  Eqs.  (5.1.2), 
(5.1.3),  (5.1.4),  (5.1.6)  and  (5.1.7)  respectively,  u and  T are  the 
expressions  (5.1.8)  and  (5.1.9)  respectively.  Also,  the  relations 
(4.1.14)  - (4.1.17)  hold  for  the  mass  transfer  problem  as  well.  The 
evaluation  of  mass  transfer  Reynolds  numbers  R^  and  R^  is  the  subject 
of  Sec.  5.3. 


5.3  Evaluation  of  Mass  Transfer  Reynolds  Numbers 

It  was  mentioned  in  Sec.  3.3  that  the  mass  transfer  Reynolds  numbers 
P.£  and  R^  have  to  be  determined  iteratively.  This  procedure  is  discussed 
briefly  in  this  section.  Referring  to  Eqs.  (3.3.6)  and  letting 


x = exp(R^Pr2)  = exp(R^)  for  Pr2 


(5.3.1) 


the  result  is 


H(x)  = Hn 


+ *n(l  - -) 
x 


I (X  - 1)  + cp(i  - ^T) 

RT  ’ T 

c (1  - ~)  + (x  - 1) 
p n i 

r x e 


rp  V 

c . T 
pi  e 


(5.3.2) 


h u2  Prl 


(5.3.3) 
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In  most  cases  of  interest  e <<  1 and  p <<  1 and  Pr^  is  typically  less 


than  10  so  that  n <<  1.  In  order  to  evaluate  R it  is  necessary  to 


solve  the  transcendental  equation  (5.3.2).  This  equation  was  solved 


using  the  Newton-Raphson  method  which  requires  specification  of  a guess 


for  x.  It  may  be  verified  that  the  function  H(x)  has  only  one  zero 


when  n <<  1 and  1/ c .T  > T /T  . The  latter  condition  holds  in  most 
pi  e w e 


practical  problems  of  interest.  Fig.  4 illustrates  the  behavior  of 


the  function  H(x)  as  a function  of  x and  it  is  seen  that  x = 1 is  the 


obvious  initial  guess.  Once  x is  determined  R is  known  from  (5.3.1) 


and  R^  is  determined  from  the  relation  R^  = epR, 


5.4  General  Solutions  of  the  Modified  Orr-Sommerf eld  Equations 


Eq.  (5.2.1)  can  be  written  in  the  form 


(^  - a^ip^"-  ~ ot^)' - “ Oj^)  “ id-jR^i^  - (u^C)  - c^) 


(5.4.1) 


= wx(C) 


(5.4.2) 


Then  Eq.  (5.4.1)  becomes 


W'l  ~ V*!  * {ai  + - c.^  } wx  = 0 
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In  order  to  reduce  the  above  equation  to  an  Airy  differential  equation 


it  is  necfcasary  to  eliminate  the  first  derivative  term.  Thus  let 


= s1(e)H1(e) 


(5.4.4) 


Substitution  of  Eq.  (5.4.4)  into  (5.4.3)  leads  to 

<2si  - Vi), . fsi  - v;  - + i“iRi<iiu) 

Hi + Hi + s2 


Eq.  (5.4.5)  indicates  that  should  be  chosen  such  that 

2S1  ~ Vl  = ° 


S1  = 6 


Rh5/2 


this  choice  of  Eq.  (5.4.5)  takes  the  form 

rv 

h'^  - Mr + ai + ia1R1^1(u  - cL)  h1  = o 


Dc  fining 


= 


— + + ia1R1(u1(4)  - c^ 


(a^u*) 


H!  = 0 


(5.4.5) 


(5.4.6) 


(5.4.7) 


(5.4.8) 


where  u^  is  a constant  and  carrying  out  the  transformation  of 
independent  variable  in  Eq . (5.4.7),  the  result  is 


(5.4.9) 
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where 


hi(ci)  * V^)) 


Eq.  (5.4.9)  is  the  Airy  differential  equation  encountered  earlier  in 
Sec.  4.3.  The  experience  gained  in  the  solution  of  the  zero  mass 
transfer  problem  suggests  the  following  solution  of  (5.4.9), 

W = c3AiU1E+)  + C4Ai(C1E“) 


HX(C)  = C3Ai(c1(5)E+)  + C4Ai(c1(C)E") 


Where  E±  is  given  by  Eq.  (4.3.12) 


From  Eq.  (5.4.4) 


RhS/2 

w2(C)  = e 11 


^C3Ai{  (C)E"t’}  + C4Ai{e1(C)E~}J 


and  hence  Eq.  (5.4.2)  assumes  the  form 


<p'l  ~ = e 


RbC/2 


C3Ai(cl(C)E+}  + C4Ai{c1(OE-} 


The  homogeneous  solution  of  the  above  equation  is 

*1H  = + C2e“al?  (5.4.] 

and  the  particular  solution  obtained  by  the  method  of  variation  of 


parameters  is 


sinh{a^(£  - t) }Ai{c3(t)E+}dt 


sinh{a^ (£  - t) }Ai{c^(t)E“}dt 
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Finally,  the  general  solution  of  Eq.  (5.2.1)  is 

c*l5  + C e'al5  + — / 

.5 


<P,  (O  = C.e 


s (v 


/2 


sinh{a^(£  - t) }Ai {4^ (t)E+}dt 


C4  [ Rh£/2 


sinh{a^(C  - t) }Ai{ 4^ (t)E- }dt 


(5.4.15) 


w.iere  C^,  C2>  C^,  and  are  arbitrary  constants  of  integration  and 
£ is  selected  such  that 


4(0  = 0 


(5.4.16) 


It  should  be  noted  that  the  same  notation  4^  is  used  in  both  zero 
mass  transfer  and  mass  transfer  cases,  however,  4^  has  different 
definitions  as  seen  from  Eqs.  (4.3.4)  and  (5.4.8). 

To  obtain  the  general  solution  of  (5.  2.3)  the  procedure  described 


above  is  repeated  to  yield  the  result 


*,(n)  = C,ea2n  + C,e"a2n  + 


— f eR&t 

a2  J 


/2 


sinh{a2  (n  - t) }Ai { 42 (t) E+}dt 


f R6c 

“l  J 6 


/2 


sinh{a2(n  - t)}Ai(42(t)E  }dt 


(5.4.17) 


where 


T + “2  + la2R2<“2<£)  " c2) 

S2(t) , “ '■'■-zR 

(a2R2U2) 


(5.4.18) 
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5.5  General  Solutions  of  Temperature  and  Concentration  Perturbation 
Equations 

(i)  Solution  of  temperature  perturbation  equation  for  liquid: 

Consider  the  homogeneous  part  of  Eq.  (5.2.2),  viz., 

e’i  - Pr-jR^  - {cx^  + ia1Pr1R1(u1(C)  - = 0 (5.5.1) 

Following  the  procedure  of  Sec.  5.4  the  first  derivative  term  is 
eliminated  from  (5.5.1)  and  the  resulting  equation  is  converted  into 
an  Airy  differential  equation.  The  latter  is  solved  as  before  to  give 
the  homogeneous  solution 

01H  = eRhPrl?/2  {C9Ai(z1(C)E'H)  + C^Ai^  (C)E“)  } (5.5.2) 


where 


zx(0 


RTPr 

— ^ — + aj  + ia^Pr^R^(u^ (O 

.2/3 


cl) 


(a1Pr1R1up" 


(5.5.3) 


Consider  now  the  particular  solution  of  Eq.  (5.2.2).  Substitution 
for  from  Eq.  (5.4.15)  into  (5.2.2)  yields 


ei  “ PriVi " ta2  + iaiPriRi(lli(c)  ' ci)ei 


= Pr  1R1'f C^l5  + C2e”al^  + 


& 


1 2 


sinh{a^(5  - t) }Ai(c^ (t)E+)dt 


ht//2  sinh{aj(^  - t)}Ai(4^(t)E  }dt( 


(5.5.4) 


where  is  constant  for  the  linear  temperature  profile  in  Eq.  (5.1.6), 
The  general  solution  of  Eq.  (5.5.4)  may  be  written 


61  = 01H  + Pr1R1T[(C1J1  + C2J2  + C3J3  + W (5.5.5) 

where  J^,  J2,  J3>  and  are  particular  integrals  obtained  by  the 
method  of  variation  of  parameters  , 


(t)y2«)  - y2(t)y  «)  E 

m e ldt 


(5.5.6) 


f y1(t)y2(0  - y2(t)yx(^)  f 

1,4  ” J WCt)  J 


Riit/2  . _ . 

e sinh{a3(E  - i)}Ai  (t)E-}dTdt 

(5.5.7) 


where 


^1,2  = eV^AiU^OE4: 


(5.5.8) 


and  the  Wronskian  of  y1  0 is  defined  by 


dy„(t)  __  dy  (t) 

w{yi(t),y2(t)}  = yx(t)  — y2(t) 


(5.5.9) 


Combining  Eqs.  (5.5.8)  and  (5.5.9)  and  simplifying,  it  may  be 


verified  that 


/1(t),y2(t)}  = eRhPr 1^  - j W{ Ai (Zl(t)E+) , Ai(Zl(t)ET)}  (5.5. 


122 


Differentiating  z^(t)  in  Eq.  (5.5.3)  w.r.t.  t 


dt 


i(alPrlRl{il)2/3 


(5.5.11) 


Also,  the  Wronskian  in  (5.5.10)  is  given  by  (Ref.  44) 

W(Ai{z.  (t)E+),  Ai  {z.  (t)E~ }]  = (5.5.12) 

l i ztr 

Introducing  the  expressions  in  Eqs.  (5.5.11)  and  (5.5.12)  into  Eq. 
(5.5.10)  one  gets 

W(t)  = ^(otlPr1R1u^)2/3eRhPrit  (5.5.13) 

It  is  new  necessary  to  evaluate  the  integrals  (5.5.6)  and  (5.5.7) 
with  the  help  cf  Eqs.  (5.5.8)  and  (5.5.13).  The  end  results  of  these 
manipulations  are  given  belcw. 

?1T>Prl£/2  ft 

J.  y= rrr  f < Ai{z  (t)E+}Ai{z  (£)E- } 

1,2  (aiPr  RGp1/3j  ( 


RhPriS/2  r 5 

2— 7 pry  fe  RhPrit/2rAi{z1(t)E+}Ai{z1(OE' 

iPriVl>  ' L 

4* 

-Ai{z1(t)E-}Al{z1(OE+}j  X 
t 

f Rh^2  ~ ~ 

le  sinh{a1(t  - t) }Ai{c^(r)E 


RhV2  _ , . „ . 

e sinh{a^(t  - x) }Ai{c^(T)E±}didt 


Finally,  the  general  solution  in  Eq.  (5.5.5)  can  be  written  in  a 
slightly  different  form  as 

RhPr i C/2 

ex(0  = e n 1 C9Ai{Zl(C)E+}  + C10Ai{Zl(^)E  } 

2irPr  R T' 

+ - ./3  ■ Vl<»  + C2J2<E)  + C3J3(t)  - V4(t) 


5. 5. 36) 


-5  RhPri 

J1  2a)  = / e(±ai  T_)t  F^tjOdt 


/ -Vrlt/2 

J3>4(0=/e 


L 

. /*  V> 

■L (t ; €)  / e 


sinh{a^(t  -i) }A1{c;^(t)E  }dxdt 


(5.5.18) 


F1(t;5)  = Ai{z1(t)E+}Ai{z1(C)E~}  - Ai{ z ± (£) E+}Ai { zx (t )E" } (5.5.19) 


r 
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and 


5x(t*)  = 0 
z}(C*)  = 0 


(5.5  20) 
(5.5.21) 


(ii)  Solution  of  temperature  perturbation  equation  in  gas-vapor: 

Following  the  procedure  described  above,  the  general  solution  of 
Eq.  (5.2.4)  is 


e2(n)  = e 


R6n/2 


Ci;LAi{z2(n)E+}  + C12Ai{z2(n)E- 


} 


2irR2T2 


* 1 / O 

(a2R2u2)i/J 


c3J5(n)  + C£J6 (n)  + c?J7(n)  + CgJg(n) 


(5.5.22) 


where 


n 

n t 

/-Rfit/2  . f R 
e F2(tjn )J  e 


e(±a2  - 2 )t  F2(t;n)dt 


(5.5.23) 


,t/2  . . _ _ 

6 sinh{a2(t  - t) }Ai(42 (x)E±}didt 


(5.5.24) 


with 


F2(t;n)  = Ai{z2(t)E+}Ai{z2(n)E-}  -Ai{ z2 (n)E+}Ai{ z2 (t)E~ } 


(5.5.25) 
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z2(t*>  = c2(n*)  = o 


It  should  be  noted  that  z^(t)  and  are  Identical  since  the 

andtl  number  for  the  gas-vapor  is  assumed  unity  (See  Eq.  (5.4.18)). 
(iii)  Solution  of  the  concentration  perturbation  equation 

Applying  the  procedure  in  (i)  to  Eq.  (5.2.5)  would  give  the 
general  solution 


BUn/2  , 

x(n)  = e c13Ai{z2(n)E+}  + C1/Ai{z9(n)E-} 


14  2 


2ttR  x j 

~ A 1 / 3 I C5J5^^  + C6J6^n^  + C7J7^  + C8J8^| 

(o2R2u2)  l J 


where  J Jg , and  Jg  are  defined  by  Eqs.  (5.5.23)  and  (5.5.24). 

The  auxiliary  equations  (5.5.25)  and  (5.5.26)  hold  in  this  case  also. 


5.6  Application  of  Boundary  Conditions 


The  boundary  conditions  (5.2.6)  - (5.2.20)  involve  derivatives  of 
^1*  ^2’  ®i»  ®2  atlC*  T^ese  differentiations  necessitate  the  use  of 
Leibniz's  rule  since  the  variables  £ and  n occur  in  the  limits  of 


integration.  Appendix  F contains  the  expressions  for  the  required 


derivatives. 
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Substituting  Eqs.  (D.l)  - (D.7),  (D.ll)  and  (D.16)  into  the 
boundary  conditions  (5.2.6)  - (5.2.10)  and  carrying  out  the  necessary 
algebraic  simplifications,  the  resulting  equations  are 

c^e-0!  C1  + o^e"1  C2  + + l^  = 0 (5.6.1) 

a,e'ai  C.  - aeai  C0  + I,C,  + I.C.  = 0 (5.6.2) 

1 1 2 3 3 4 4 

a2ea2  C5  + a2e_a2  C(>  + 1 5C?  + I6Cg  = 0 (5.6.3) 

a2ea2  C5  - a2e"a2  Cg  + I?C7  + I?Cg  = 0 (5.6.4) 


ualCl  ~ ualC2  + milC3  + mi2CU  ~ £a2C5  + ea2C6  ' eI15C7  " eI16C8 

R R. 

= i (uu|  - eu2)  + ^(^R a^u  — ) (5.6.5) 

2ua^C1  + 2ua2C2  + | uAi{c1(0)E+}  - 20^  J C3 

+ [ uAi{<;1(0)E_}  -2a1I10  J C4  - 2TTe2C5  - 2iJe2C6 

+ j--ye2Ai{c2(0)E+}  + 2Ue2a2I13  | C?  + j-Ue2AiU2(0)E"}  + 2u£2a2I14  [ Cg 

= 0 (5.6.6) 
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- V“i)l  ci 

2 | Rv,  A / \)  2TIq11 

+ -f|^^(ai  + (lii(0)-ci)aij|  -wdc2 
•“  ' 

+ — I - Ai{;1(0)E+}  + Ai,k,(0)E+k;E+  . 

P a2  [ 2 Ri  1 Rx  1 

+ £■(——  + ■—  • I +(—  £ (u  (0)  - c ) + 4^—1  I C 

P Ri  9 |ai  p 1 1 £TrR2j  11  3 

+ - £■  V'  " TIT  A1(C1(0)E-}  + i-  Ai'{C1(0)E+)c;E-  • 

p a2  Z Rx  1 J-  1 


F 


a1a2ilC1  + axa2uC2  - - c^uI^C^ 


ala2^^5  ~ ala2p^6  + alp^'13^7  + alp^'14^8 
“ ia^2[  G(G1(0)  - ci)  - p"(u2(0)  - c2)] 

J^R2 C1  " X(0))+  P2{z2I3l  ~ 2~  I295]  C5 

+ jl^ (1  - X(0))+  P2^2J'3?  - 2~  ^30^  C6 

+ E”  “2(1  **  ^(0))I13+  P2{z2I47  “ 2_I45}]  C7 

+ ^2<'1  ” ^0))I1A+  P2{z2I48  “ ~I46}]  C8 

+ 2^Ai{z2(0)E+}  + z2Ai’{z2(G)E+}E+j  C13 

+ ^Ai{z2(0)E"}  + z2Ai'{z2(0)E"}E_J 


= ia^R2 (1  _ X(0))(u2(0)  - c2) 


(5.6.8) 


(5.6.9) 


where 


2wR2x 

P2  ~ , n >■  vl/3 
(a2R2u2) 


Q1I17C1  + Q1I18C2  + Q1X33C3  + Q1I34C4 


+ Ai{z1(-1)E+}C9  + Ai{z1(-1)E"}C10  = 0 (5.6.10) 
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*1  = 


2-irPr^T^ 

" . 1/3 

(a1Pr1R1u1) 


where 


Q2I21C5  + Q2I22C6  + Q2I37C7  + Q2I38C8 


+ Ai( z2 (1)E+}C11  + Ai{z2(l)E"}C12  = 0 


2ttR2T2 


(5.6.11) 


2 / B * W3 

'a2^2u2 ^ 


QlTI25Cl  + QlTI26C2  + Q1TI41C3  + Q1TI42C4 

-Q2I29C5  " Q2I30Co  " Q2I45C7  " Q2r46C8 
+ TAi{z1(0)E+}C9  + TAi{z1(0)E"}C10 
- A1{z2(0)E+}C11  - Ai{z2(0)E~}C12  = eT2  - TTj. 


TQ 

, VrL  1 

Cl  + 

TQ 

„ * Vri , I 

~ 

eZ 

^^27  + 2 :25  J 

“ etr  ' 
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[xW  A^V0^]  C13  + (jW  Ai{z2(0)E"}  j = 0 


- Ai{z2(0)E~}  C12 


where 


y2  = -fro) 


(t2(0)} 


The  definitions  of  integrals  1^  - I^g  are  listed  in  Appendix  E. 

As  in  the  case  of  zero  mass  transfer  the  boundary  conditions 
(5.6.1)  - (5.6.7)  and  (5.6.9)  - (5.6.15)  form  a system  of  linear 
algebraic  equations  of  the  type 

[A(c1)]{c)  = {V(C;L) } (5.6 

where (A)is  a 14  x 14  coefficient  matrix  of  the  left  hand  sides  and 
(V(c^) } is  a 14  x 1 column  vector  of  the  right  hand  sides.  The 
remaining  equation  represents  the  characteristic  function  for  the 
mass  transfer  problem  and  can  be  written  in  the  form 

G jc^C^^lJ  = a2u  - I9C3  - I10c4j 

“V  j^a2C5  + a2C6  " I13C7  " 1 l4C8^j 
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-ia^a^  u 


(ux(0) 


cx)  - P(u2(0) 


c2) 


)j  - 0 


(5.6.17) 


or  in  a condensed  notation 


^c1; {C(c1>  }J  - (CCc^  }T{U(c1) } -ia^a2^  u(u1(0)  - c - p(u2(0)  - c2)J 

(5.6.18) 


where 


and 


-1 


(C^)}  = (A(Cl)  ) {V(cT)  } 


(5.6.19) 


{U(c1)}T=  ^aia2u  o‘1ot2u  -a2uI9  -a2uI10 

-a1a2p"  aipIi3  alpI14^  (5.6.20) 

In  Eqs.  (5.6.17)  and  (5.6.18),  c2  is  again  given  by  Eq.  (4.1.15). 

As  in  Chapter  IV  G is  the  characteristic  function  and  is  the 
eigenvalue.  The  function  G can  be  expressed  as 

G(aL,  c,v,  k,  P,  c R2,  W.F,  E,  R^,  R^ , Pr^  A,  R;  c±)  = 0 (5.6.21) 

The  stability  problem  is  therefore  reduced  to  locating  the  zeros  of 
an  analytic  function  G in  the  complex  c^  plane.  The  method  of  zero 
finding  is  described  in  the  next  section. 

5 . 7 Outline  of  the  Eigenvalue  Iteration  Procedure 


The  zeros  of  the  characteristic  function  G are  determined 
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numberically  using  methods  similar  to  those  of  Sec.  4.5.  The 
important  steps  are  listed  below. 

(i)  The  initial  guess  for  the  eigenvalue  is  obtained  from  the 
solution  of  the  zero  mass  transfer  problem.  This  restricts  the  inves- 
tigation to  the  effects  of  mass  transfer  on  particular  modes  and  the 
question  whether  mass  transfer  itself  introduces  any  stability  modes 
is  left  unanswered.  The  latter  statement  will  become  more  meaningful 
in  conjunction  with  the  developments  of  Chapter  VI. 

(ii)  For  given  parameter  values  in  Eq.  (5.b.21)  and  the  guess  value 

c^,  evaluate  the  integrals  1^  through  I^g  using  the  integration  procedures 
described  in  Appendices  E and  G.  -* 

(iii)  Generate  the  coefficient  matrix  (A)  and  obtain  its  inverse.  The 
IBM  routine  MINV  was  used  for  this  purpose  with  minor  modifications. 

Also  generate  the  column  vector  {V}. 

(iv)  Determine  the  constants  of  integration  {C}  by  computing  the 
product  (A]  ^{V}. 

(v)  Calculate  G in  Eq.  (5.6.17).  Usually  this  equation  will  not  be 
satisfied  with  the  first  guess  for  c . 

(vi)  Obtain  an  improved  approximation  for  c^  by  employing  the  Newton- 
Raphson  iteration  technique  described  in  Appendix  I. 

(vii)  Compare  successive  values  of  c^  for  convergence  within  a 
prescribed  tolerance  on  real  and  imaginary  parts.  Repeat  steps  (ii)  - 
(vi)  until  desired  convergence  is  reached. 


6 . 1 Description  of  Data 

22 

The  experimental  data  of  Craik  was  chosen  for  numerical  compu- 
tations. The  present  theoretical  model  does  not  match  the  experimental 

% 

conditions  perfectly  and  hence  some  of  the  parameters  in  Craik' s inves- 
tigation were  suitably  ' corrected' . For  instance,  the  experiments  were 
conducted  in  a closed  rectangular  channel  1 inch  high  and  6 inches  wide. 
The  air  flew  was  provided  by  a large  fan  which  drew  air  through  the 
apparatus.  Water  was  introduced  into  the  channel  at  the  entry  section 
and  formed  a film  on  the  bottom  made  of  plate  glass.  Consequently  the 
velocity  profile  in  air  was  parabolic  but  the  liquid  velocity  profile 
was  very  nearly  linear.  It  was  assumed  in  the  present  analysis  that 
the  velocity  profile  in  the  gas  is  linear  and  has  a thickness  6.  There- 
fore, in  order  to  adapt  Craik' s data  to  the  present  model,  <5  was  chosen 
to  be  half  the  difference  between  channel  height  and  liquid  film  thick- 
ness. The  velocity  profile  between  the  interface  and  y =8  was  assumed 
linear. 

Another  important  correction  to  Craik' s data  was  made  in  regard  to 
the  value  of  gas  viscosity.  The  air  flow  in  his  experiments  was  turbu- 
lent and  hence  an  augmented  value  of  laminar  viscosity  must  be  used. 

This  is  accomplished  as  follows.  Craik  measures  the  steady-state  inter- 
face velocity  u^  and  determines  the  film  thickness  h by  measuring  the 
» volumetric  flow  rate  of  liquid.  This  permits  one  to  compute  the  inter- 

facial shear  cn  the  liquid  side  using  the  expression 
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Vif 


(6.1.1) 


Since  this  value  of  shear  stress  must  equal  the  interfacial  shear  on 
the  gas  side,  it  follows  that 

y2Ue 

tx  » t2  = -j-S-  (6.1.2) 

and  therefore, 

y2*h~Ml  (6-1-3) 

e 

Since  the  gas  Prandtl  number  is  assumed  unity,  the  thermal  con- 
ductivity k.£  of  the  gas  for  turbulent  flow  is  calculated  using  the 
relation 


k2  = 


k2 iam^k 


(6.1.4) 


where  the  factor  is  given  by 


F 


k 


^2 

^2  S.am 


(6.1.5) 


In  Eqs.  (6.1.4)  and  (6.1.5)  anc*  ^enote  molecular  thermal 

conductivity  and  viscosity  respectively. 

The  list  of  various  dimensional  and  non-dimensional  parameters 
used  in  numerical  computations  is  contained  in  Table  II  at  the  end  of 
this  chapter.  The  data  in  this  Table  corresponds  to  both  air  and 


water  at  room  temperature. 
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6,2  Root  Location  Procedure  for  Zero  Mass  Transfer  Problem 


One  of  the  eigenvalues  (or  modes)  can  be  immediately  calculated 
from  the  long  wavelength  disturbance  solution  of  Sec.  4.2.  Thus  sub- 
stituting for  e and  p from  Table  II  into  the  approximate  expression 
(4.2.19)  and  the  exact  expression  (4.2.18)  the  results  are 


c,  = 1.086  + 01 

lapprox 


c.  _ = 1.064  + Oi 

lexact 


for  <<  1 and  <<  1 


Starting  with  the  guess  c^  = 1.086  and  choosing  = 0.001  (hence 
a2  = ai/e  = 0*023)  the  Newton- Raphson  procedure  yields  the  result 

c1  = 1.06438  - 0.00774  for  = 0.001 

Once  c^  is  known  for  a given  a^,  it  is  a simple  matter  to  trace  this 
mode  by  varying  gradually. 

Identification  of  other  stability  modes  is  somewhat  complicated 
and  tedious.  The  reader  is  referred  to  Sec.  4.6  in  this  connection. 

An  illustrative  example  of  approximate  location  of  roots  of  the  charac- 
teristic equation  (4.4.11)  is  given  in  Figs.  5 and  6.  Suppose  it  is 
desired  to  find  the  eigenvalues  (i.e.  roots  of  the  characteristic 
equation)  in  the  interval  1 < c^r_S  2 and  -1  _<  c^  _<  0.  Then  the 
function  G in  Eq.  (4.4.16)  is  calculated  at  c^  = 0,  -0.1,  -0.2,  ..., 
-1.0  for  each  c^r  = 1.0,  1.1,  ...,  2.0.  The  next  step  is  to  plot 
Re(G)  and  Im(G)  against  c^with  c^  as  a parameter.  The  results  are 


; 
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shown  in  Figs.  5 and  6 for  the  case  a ^ = 0.05.  This  choice  of  was 
dictated  by  the  fact  that  the  system  is  expected  to  be  well-behaved  for 
long  wavelength  disturbances.  It  is  seen  that  both  Re(G)  and  Im(G)  go 
to  zero  around  c^  = 1.32  - 0.4i  and  c^  = 1.84  - 0.5i.  With  these 
guess  values  the  Newton-Raphson  iteration  gives  the  exact  results 

c1  = 1.318  - 0. 405i 

and  for  ai  = °-05 

c1  = 1.842  - 0. 505i 

This  process  of  root  location  was  carried  out  in  the  interval 
0 _<  c^^  _<  4 and  -2  c^  < 2.  The  following  spectrum  of  eigenvalues, 
ordered  according  to  real  part,  was  obtained: 

For  = 0.05, 

cn  = 0.273  - 0. 197i 
c12  = 1.318  - 0. 405i 

C13  = 1,411  + 1’053i 

c. . = 1.842  - 0.5051 
14 

c15  = 3.263  - 1 . 220i 

Note  that  only  c^  has  a critical  point  inside  the  liquid  (i.e., 
c^r  < 1.)  A random  search  for  an  eigenvalue  with  c^r  > 4 led  to  the 
root 


17.48  + 0.5401 
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It  should  be  emphasized  that  this  eigenvalue  was  obtained  by  supplying 
random  guesses  (with  > 4)  and  was  not  obtained  through  the 
systematic  search  procedure  mentioned  earlier. 


6.3  Amplif ication  and  Phase  Velocity  Curves  for  Zero  Mass  Transfer  Case 


It  was  stated  in  Chapter  I that  the  stability  or  instability  of 
the  interface  depends  upon  whether  an  infinitesimal  disturbance  of  a 


given  wavelength  grows  or  decays  with  time.  As  pointed  out  in  Sec.  2.7 
positive  (or  c^)  corresponds  to  an  unstable  interface  and  negative 
corresponds  to  a stable  interface.  With  these  comments  in  mind  one 
seeks  to  know  hew  (or  more  correctly  a^c.^)  varies  with  the  distur- 
bance wave  number  k(=  2tt/A).  This  requires  that  the  modes  obtained  in 
Sec.  6.2  for  = 0.05  be  traced  as  changes  continuously.  Each  of 
the  six  modes  mentioned  in  the  previous  section  was  traced  by  varying 
very  gradually.  Much  care  needs  to  be  exercised  during  this  process 
since  a sufficiently  large  change  in  may  result  in  switching  to  a 
different  mode.  This  exercise  was  carefully  done  and  the  results  are 
presented  in  Figs.  6 through  11  as  amplification  and  phase  velocity 


curves.  The  plots  of  phase  velocity  against  wave  number  have  been 
included  for  the  sake  of  completeness  and  also  to  facilitate  comparisons 
with  the  well-known  water  wave  phenomena  such  as  gravity-surface  tension 
waves  and  Kelvin-Helmholtz  waves.  The  curves  of  amplification  have  been 
presented  in  the  form  a^c.^  vs  and  the  phase  velocities  are  plotted 
in  the  form  a.^c^r  vs  It  may  he  recalled  that  a^c^  and  aicqr  are 


A very  important  point  needs  to  be  brought  to  the  attention  of 
the  reader.  Only  those  modes  discovered  for  = 0.05  in  the  previous 
section  have  been  traced  as  varies.  It  is  conceivable  that  more 
stability  modes  ’creep  in'  as  increases.  In  fact,  some  evidence 
was  obtained  during  the  numerical  investigation  which  showed  that 
this  is  true.  It  should  soon  become  clear  that  the  six  modes  in  the 
present  investigations  have  rather  distinct  characteristics.  These 
are  discussed  below. 

(i)  Figs.  6a  and  6b  represent  amplification  and  phase  velocity 
characteristics  for  the  eigenvalue  c^..  This  eigenvalue  is  of  some 
interest  because  it  is  the  only  one  with  a critical  point  inside  the 
liquid.  Fig.  6a  shows  that  the  imaginary  part  of  c is  always  less 
than  zero  and  hence  this  mode  is  stable  for  all  a^.  The  curious  dip 
around  = 0.65  is  not  explained.  The  phase  velocity  plot  of  Fig.  6b 
also  exhibits  sharp  changes  at  this  value  of  a^.  The  phase  velocity  is 
generally  increasing  with  a^. 

(ii)  The  eigenvalue  c^  is  interesting  because  it  was  the  only  one 
predicted  analytically.  It  may  be  recalled  from  Sec.  4.2  that  this 
mode  (at  least  for  small  values  of  a^)  has  a physical  interpretation 
that  it  depends  only  on  the  thickness  ratio  c and  viscosity  ratio  y", 
and  it  is  independent  of  Reynolds,  Froude  and  Weber  numbers.  The 
amplification  curve  in  Fig.  7a  shows  that  this  eigenvalue  is  also 
stable  for  all  ot^  and  displays  a small  dip  around  = 0.6.  The 
foregoing  observations  suggest  that  this  eigenvalue  can  be  associated 
with  the  Tolimien-Schlichting  mode  of  stability  (Sec.  1.2i)>  The 
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phase  velocity  (Fig.  7b)  increases  continuously  with  but  at  a 

♦ 

much  faster  rate  than  the  c^  mode. 

(iii)  The  eigenvalue  c^  was  found  to  be  the  only  unstable  mode  at 

» al  = ^-05  and  deserves  attention  for  this  reason.  The  variation  of 

amplification  rate  (Fig.  8a)  shews  that  this  mode  is  unstable  for  all 
(except  at  = 0 where  it  is  neutrally  stable).  It  is  seen  that 
the  rate  of  amplification  increases  rapidly  beyond  - 0.5.  The 
phase  velocity  plot  of  Fig.  8b  also  displays  a sharp  change  around 
this  value  of  a^.  The  phase  velocity  increases  with  in  this  case 
also. 

(iv)  Fig.  9a  represents  amplification  curve  for  the  eigenvalue  c.^. 

It  exhibits  a unique  characteristic  in  that  this  mode  is  stable  at 
small  values  of  and  unstable  at  large  a ^ . In  other  words.  , there 
are  two  distinct  regions,  one  stable  and  the  other  unstable,  separated 
by  a neutral  stability  point.  With  the  exception  of  c^  (which  will 

be  discussed  shortly)  none  of  the  other  modes  demonstrate  this  behavior. 

It  appears  that  this,  mode  is  the  same  as  the  one  obtained  by  Bordner, 

37  30 

et.  al  using  the  data  of  Cohen  and  Hanratty  . It  is  interesting  to 

compare  the  phase  velocity  for  this  eigenvalue  with  the  speed  of 

propagation  of  surface  tension  - gravity  waves  and  Kelvin-Helmholtz 

waves  (Fig.  9b).  It  is  observed  that  long  wavelength  disturbances 

(i.e.,  small  a^)  of  this  mode  propagate  with  nearly  the  same  speed  as 

Kelvin-Helmholtz  waves.  The  mode  under  consideration  can  be  asso- 

9 

ciated  with  the  class  C (or  Kelvin-Helmholtz)  instability  of  Benjamin 
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and  Landahl^.  This  eigenvalue,  therefore,  is  called  the  modified 
Kelvin-Helmholtz  mode  in  the  present  investigation.  The  expressions 
used  for  calculating  the  speeds  of  surface  tension-gravity  waves  and 
Kelvin-Helmholtz  waves,  in  dimensionless  form,  are  given  below: 


c2 

o 


1 - p + a2W2F2 
coth(a^)  + TX  cothCa^) 


surface  tension- gravity  wave 


(6.3.1) 


coth(a^)  + 1=  coth(a2) 
coth(a^)  + p"  coth(a2) 


^ j 1 - p + a2W2F2  | 

a f2  | coth(a^)  + p"  coth(a2)j 


1/2 


p (1  - =)2  coth(a1)coth(a2) 
--  ” 

(cothCa^)  + p coth(a2>} 


Kelvin-Helmholtz  wave  (6.3.2) 

These  equations  hold  for  two  fluids  bounded  between  two  walls  at 
y = -h  and  y = 6 respectively. 

Finally,  it  is  observed  from  Fig.  9a  that  instability  sets  in 
when  is  above  0.145.  Corresponding  to  this  value  of  a^,  = 1.99 

from  Fig.  9b.  The  dimensional  value  of  c^  is  then  13.5  cm/ sec  and 
compares  favorably  with  Craik's  experimental  value  of  11.9  cm/sec  at 
the  onset  of  instability. 

(v)  The  amplification  and  phase  velocity  curves  corresponding  to  the 
eigenvalue  c^,.  are  plotted  in  Figs.  10a  and  10b  respectively.  The 
phase  velocity  increases  with  similar  to  the  other  modes  and  complete 
stability  prevails  for  all  as  shown  by  the  amplification  plot.  Thus 
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this  mode  does  not  display  any  peculiar  characteristics. 

(vi)  The  eigenvalue  c^  is  a representative  fast  moving  disturbance 
and  therefore  merits  some  consideration.  Fig.  11a  shews  that  the  ampli- 
fication curve  oscillates  rapidly  with  and  is  somewhat  irregular. 

Thus  there  are  several  regions  of  stability  and  instability  separated 
by  neutral  stability  points.  The  phase  velocity  curve  (Fig.  lib)  also 
exhibits  irregular  behavior.  In  fact,  the  phase  velocity  decreases 
with  a,  over  a small  range. 

The  above  description  of  six  different  stability  modes  shews  that 
amongst  the  slow  moving  disturbances  the  modified  Kelvin-Helmholtz 
mode  is  the  most  interesting.  Therefore,  further  attention  is  concen- 
trated on  this  particular  mode. 

6 . 4 Effect  of  Neglecting  Instability  in  Gas 

It  was  me  itioned  in  Secs.  1.2  and  2.1(2)  that  the  validity  of  the 
frequently-made  assumption  of  neglecting  the  phase  speed  in  gas  distur- 
bance equations  is  doubtful.  To  examine  this  question  the  zero  mass 
transfer  problem  was  solved  putting  c^  = 0 in  gas  disturbance  equations 
(4.1.2),  (4.1.9)  and  (4.1.11).  This  amounts  to  assuming  the  disturbance 
in  Eq.  (2.6.7)  to  be  of  the  form 

q21(x»y»t)  = q2(y)elkx  (6.4.1) 

This  exercise  required  only  minor  modifications  in  the  Newton-Raphson 
procedure  and  in  the  computer  program.  As  pointed  out  in  the  previous 
section,  only  the  modified  Kelvin-Helmholtz  mode  was  considered.  The 
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results  of  these  computations  are  shown  in  Figs.  12a  and  12b,  again 
in  the  form  of  amplificati jii  and  phase  velocity  curves.  The  amplifi- 
cation plot  for  c^  0 is  a portion  of  Fig.  9a.  A comparison  of  the 
curves  for  C£  1 0 and  c.^  = 0 reveals  the  following  significant  results. 

(i)  The  assumption  of  neglecting  the  instability  in  gas  has  no 
effect  at  low  wave  numbers  (i.e.  for  disturbances  with  long  wavelengths) 
and  affects  the  neutrally  stable  wave  number  only  slightly. 

(ii)  Beyond  the  neutrally  stable  wave  number  the  above  assumption 
results  in  underestimation  of  the  amplification  rate. 

(iii)  When  the  wave  number  is  sufficiently  large  (e.g.  a > 0.4  in 
Fig.  12a)  the  c.^  = 0 assumption  predicts  stability  when  the  interface 
is  actually  unstable. 

Fig.  12b  shows  that  there  is  no  appreciable  difference  between 
phase  velocity  curves  for  the  cases  c^  ¥ 0 and  c ^ = 0.  Another  impor- 
tant aspect  of  the  assumption  under  question  needs  to  be  investigated. 

21 

It  was  pointed  out  in  Secs.  1.2  and  2.1(2)  that  Benjamin  established 
the  criterion  which  allows  one  to  make  the  rigid  wavy  wall  (or  C£  = 0) 
assumption.  This  suggests  calculation  of  Benjamin's  parameter  in  Eq. 
(2.1.1)  for  different  values  of  a^.  This  parameter  written  for  the 
gas  (i.e.  fluid  '2')  is 


B 

P 


m2C2r 


« 1 


(6.4.2) 


Fig.  13  shows  a plot  of  against  . It  is  seen  that  <<  1 
holds  only  for  very  small  values  of  (typically  < 0.1).  When 
> 0.2  Benjamin's  parameter  can  no  longer  be  considered  small  com- 
pared to  unity.  This  discussion  explains  why  the  amplification  curves 
for  c.^  4 0 and  c.^  ~ 0 display  characteristically  different  behavior 
when  > 0.1. 

6.5  Effects  of  Interface  Mass  Transfer 

A typical  example  of  the  effects  of  interface  evaporation  on  the 
modified  Kelvin- Helmholtz  mode  is  shown  in  Figs.  14a  and  14b.  Craik's 
data  in  Table  II  have  been  used  in  these  computations.  It  should  be 
noted  that  the  evaporative  mass  transfer  occurs  at  room  temperature 
(68  F)  and  consequently  the  mass  transfer  Reynolds  numbers  and  R^ 
are  very  small.  The  amplification  and  phase  velocity  curves  with  and 
without  mass  transfer  are  included  in  Figs.  14a  and  b.  It  is  observed 
that  both  amplification  curves  coincide  in  the  stable  range  and  the 
neutrally  stable  wave  number  is  very  slightly  affected.  When  is 
beyond  the  neutrally  stable  value,  however,  interface  mass  transfer 
results  in  an  increase  in  amplification  rate.  It  can  therefore  be 
concluded  that  the  effect  of  mass  transfer  is  destabilizing  for  this 


particular  stability  mode. 
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A comparison  of  phase  velocity  curves  in  Fig.  14b  shows  that 
mass  transfer  leads  to  a small  increase  in  the  phase  velocity.  In 
both  Figs.  14a  & b,  the  mass  transfer  effects  seem  to  become  more 
significant  as  increases.  This  is  understandable  since  a highly 
rippled  interface  (ot^  large)  causes  increased  mass  transfer  pertur- 
bations. 

The  mass  transfer  curves  are  terminated  at  = 0.25  because 
the  Newton-Raphson  procedure  failed  to  converge  beyond  this  point. 
Further  solutions  were  attempted  using  a conjugate  gradient  method 
(a  variation  of  the  method  of  steepest  descent) . This  method  is 
extremely  slow  because  it  calls  the  mass  transfer  program  a large 
number  of  times.  The  workability  of  this  method  was  tested  for 
< 0.25  during  the  last  stages  of  the  present  investigation. 
Therefore,  eigenvalue  solutions  for  > 0.25  could  not  be  included 
in  the  present  report. 

6.6  Suggestions  for  Future  Investigations 

The  experience  gained  in  the  present  investigation  suggests  that 

the  following  areas  of  the  stability  problem  need  further  study. 

(i)  The  different  modes  of  stability  should  be  classified  according 

9 1 

to  their  physical  interpretations.  The  works  of  Benjamin  and  Landahl 
should  serve  as  models  in  this  quest.  It  appears  that  a stability  mode 
associated  with  a particular  physical  mechanism  could  be  'singled  out' 


14/ 


from  the  rest  by  employing  some  approximate  technique.  An  insight 
into  the  physical  phenomena  is  necessary  to  understand  which  modes 
occur  in  practice. 

(ii)  A combination  of  physical  and  mathematical  reasoning  is  re- 

quired  to  understand  the  structure  of  the  entire  eigenvalue  spectrum. 

13 

Recent  investigations  of  Mack  show  that  the  eigenvalue  spectrum  of 
a laminar  boundary  layer,  with  a discontinuous  first  derivative  in 
velocity,  is  infinite. 

(iii)  The  work  on  the  liquid  film  stability  problem  must  eventually 

2 

include  the  effects  of  mean  velocity  profile  curvature.  Miles'  work 
shows  that  velocity  profile  curvature  plays  an  important  role  in  the 
energy  transfer  mechanism.  Since  the  Orr-Sommerfeld  equation  cannot 
be  solved  exactly  for  a general  velocity  profile,  a finite  difference 
approach  will  have  to  be  adopted.  Therefore,  one  may  consider  solving 
the  present  linear  velocity  profile  problem  using  finite  differences  as 
a first  step.  This  exercise  will  help  in  gaining  familiarity  with  the 
problems  of  eigenvalue  location. 

(iv)  The  stability  problem  with  mass  transfer  needs  to  be  studied  in 
greater  detail,  even  for  an  incompressible  air  flow.  For  instance,  the 
effect  of  mass  transfer  on  different  stability  modes  should  be  investi- 
gated. Further  numerical  investigations  should  be  carried  out  with  the 
present  model  for  higher  rates  of  mass  transfer.  An  important  aspect 
of  this  problem,  not  covered  in  the  present  work,  is  whether  the  physi- 
cal process  of  mass  transfer  itself  introduces  any  modes  of  instability. 
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(v)  Finally,  one  wishes  to  solve  the  problem  of  compressible  boundary 
layer  over  a liquid  film.  The  processes  of  heat  and  mass  transfer  will 
undoubtedly  be  extremely  important,  especially  when  the  air  flow  is 
* supersonic.  The  latter  case  is  interesting  because  multiple  stability 

loops  are  known  to  exist  when  the  mean  flow  in  the  boundary  layer 
becomes  supersonic. 


TABLE  II 

LIQUID  PROPERTIES  AT  528  DEG  R 

Liquid  = Water 
Molecular  Weight  = 18.0 

Latent  Heat  of  Vaporization  = 971.65  Btu/lbm. 

Coefficient  of  Viscosity  = 2.107  x 10  ^ lbm-sec/ft2. 

Density  = 1.937  slug/ft3. 

Specific  Heat  = 1.00  Btu/lbm-deg  R. 

Thermal  Conductivity  = 9.611  x 10  5 Btu/f t-sec-deg  R. 

Surface  Tension  = 4.926  x 10  3 lbf/ft. 

_ 3 

Liquid  Layer  Thickness  • 1.755  x 10  ft. 

Wall  Temperature  = 528.0  deg  R. 

GAS  PROPERTIES  AT  528  DEG  R 
Gas  = Air 

Coefficient  of  Viscosity  = 5.549  x 10  6 lbf/ft2. 

Density  = 2.340  x 10  3 slug/ft3. 

Specific  Heat  = 0.24  Btu/lbm-deg  R. 

Thermal  Conductivity  = 6.100  x 10  5 Btu/f t-sec-deg  R. 

Velocity  at  Edge  of  Boundary  Layer  = 19.66  ft/sec. 

Static  Pressure  = 2.116  x 103  lbf/ft2. 

Temperature  at  Edge  of  Boundary  Layer  = 528.0  deg  R. 

-2 

Boundary  Layer  Thickness  = 4.079  x 10  ft. 

CHARACTERISTIC  NON-DIMENSIONAL  PARAMETERS 

_ 2 

Liquid  Layer /Boundary  Layer  Thickness:  e = 4.302  x 10 

Gas  Viscosity /Liquid  Viscosity:  "y  = 0.263 

_ - 3 

Gas  Density/Liquid  Density:  p = 1.208  x 10 

Gas  Specific  Heat/Liquid  Specific  Heat:  c^  = 0.240 

Gas  Thermal  Conductivity/Liquid  Thermal  Conductivity:  k = 0.635 
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Liquid  Reynolds  Number:  = 35.54 

Gas  Reynolds  Number:  R^  = 338.2 

Mass  Transfer  Reynolds  Number  for  Liquid:  R^  = 2.95  x 10 

Mass  Transfer  Reynolds  Number  for  Gas:  R^  = 2.60  x 10  2 

Liquid  Prandtl  Number:  Pr^  = 10.04 

Gas  Prandtl  Number:  P^  = 1.0 

Liquid  Weber  Number:  W = 5.46 

Liquid  Froude  Number:  F = 0.93 

Gas  Euler  Number:  E = 2.34  x 103 


152 


The  following  general  conclusions  can  be  drawn  from  the  present 
analysis  of  liquid  film  stability. 

(i)  The  stability  of  an  interface  between  two  fluids  is  charac- 
terized by  the  existence  of  several  modes  of  stability.  These  modes 
may  be  completely  stable,  unstable  or  may  change  from  stable  to  un- 
stable (or  vice  versa)  as  the  disturbance  wave  number  is  varied.  Two 
modes  were  distinguished  phvsically  in  this  investigation,  one  asso- 
ciated with  the  Tollmien-Schlichting  instability  and  the  other  with 
Kelvin-Helmholtz  instability.  It  should  be  pointed  out  that  the 
foregoing  observations  hold  for  the  case  of  linear  velocity  profiles 
in  both  the  gas  and  the  liquid.  The  assumption  of  linear  profiles  is 
introduced  in  order  to  isolate  the  effects  of  mass  transfer  on  stability. 
Further  research  is  recommended  to  study  the  eigenvalue  spectrum  for 
curved  velocity  profiles. 

(ii)  The  customary  assumption  of  neglecting  instabilities  in  the 

incompressible  gas  motion  is  valid  only  for  very  small  values  of  the 
disturbance  wave  number  (a^  <<  1) . When  the  disturbance  wave  number 
is  moderate,  i.e.,  = 0(1),  such  an  assumption  not  only  leads  to  a 

gross  underestimation  of  the  amplification  rate  (for  the  modified 
Kelvin-Helmholtz  mode)  but  can  even  predict  incorrectly  a stable 
interface.  Again,  these  conclusions  apply  to  linear  velocity  profiles 
in  both  fluids  and  need  to  be  extended  to  include  the  effects  of 


velocity  profile  curvature. 
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(lii)  Limited  computations,  for  very  small  mass  transfer  rates, 
indicate  that  interface  evaporation  has  a negligible  effect  on  the 
modified  Kelvin-Helmholtz  stability  mode  for  very  small  disturbance 
wave  numbers  and  has  a destabilizing  effect  when  the  wave  number  is 
moderate. 
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APPENDIX  A 

ENERGY  EQUATION  FOR  LIQUID 

The  equation  of  state  for  a liquid  treated  as  a pure  substance 
is  of  the  functional  form 


hl  = hi(^pi’Tl) 


(A.  1) 


where  the  notation  of  Chapter  II  is  preserved.  A small  change  in 
the  enthalpy  h^  for  a closed  system  can  be  expressed  as 


dh 


3 

1 9p 


»h. 


dp  + 


1 3T, 


dT, 


(A. 2) 


Modifying  Eq.  (A. 2)  for  a differential  control  volume  moving  with  a 
fluid  particle  the  result  (stated  without  proof)  is 


Dh,  3h, 


Dt  3p\ 


Dp1  31k 
Dt~  + 3T^ 
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(A.  3) 


For  a pure  substance 
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and 


£7  (1  - W 


(A.  5) 


where  6,  Is  the  coefficient  of  volume  expansion  and  c , is  the  specific 
1 P1 

heat  at  constant  pressure.  Substituting  Eqs.  (A. 4)  and  (A. 5)  into  Eq. 
(A. 3),  using  the  definitions  of  substantial  derivatives  and  combining 
with  Eq.  (2.3.4)',  the  result  is 


DTj 

picPi  dT 


_ DPi 

'6lTl^T+klV2Tl 


(A. 6) 


where  viscous  dissipation  is  neglected.  The  discussion  in  Chapter  I 
indicates  that  the  pressure  gradient  in  the  present  problem  should  be 
small.  Also,  the  term  3^T^  is  small  (e.g.  for  water  3^  = 304  x 10  6/ 
deg  C and  for  in  the  range  of  100  deg  C,  3^T^  ~ 0.03)  and  there- 
fore the  first  term  on  the  right  hand  side  of  Eq.  (A. 6)  can  be  ne- 
glected. This  equation  then  reduces  to 


DT 

picPi  dT  = ki72Ti 


(A. 7) 


which  is  a compact  form  of  Eq.  (2.3.4). 
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APPENDIX  B 


DERIVATION  OF  SHEAR  AND  NORMAL  STRESS  EQUATIONS 


Consider  equilibrium  of  a triangular  element  of  fluid  at  the 
interface  (Fig.  B.l).  Resolving  all  the  forces  normal  and  tangen- 
tial to  As  and  summing  up,  the  following  expressions  are  obtained: 

a = o sin2*  + o cos2*  - 2t  sin*cos4>  (B.l) 

xx  yy  xy 

t ■ y(o  - o )sin2if>  + t (cos2<f>  - sin2<J>)  (B.2) 

» yy  xx  xy 


For  an  incompressible  fluid  the  stress  tensor  is  given  by 


_ 3u 

°xx  = "P  + 2p  77 


„ , 0 8v 

o = -p  + 2y  -x — 

yy  K 3y 


(b.  3) 


T = V 
xy 


8u  8v 
[8y  8xJ 


From  the  geometry  of  Fig.  B.l,  tan<}>  = n , so  that 


sin$ 


(1  + nj) 


1/2 


and 


(B.  4) 


« + n5)1/2 


COS* 


Introduction  of  Eqs.  (B.3)  and  (B.4)  into  Eqs.  (B.l)  and  (B.2)  yields, 
after  some  simplifications. 


These  expressions  are  utilized  in  Chapters  II  and  III. 


167 


APPENDIX  C 


DERIVATIVES  OF  GENERAL  SOLUTIONS  WITHOUT  MASS  TRANSFER 

Derivatives  of  (£)  and  ^(n)  in  Eqs.  (4.3.17)  and  (4.3.19) 
w.r.t  £ and  n,  obtained  using  Leibniz's  rule  are 

*K>-c 


+ C 


4 

^ cosh{a^(£,  - t)  }Ai(c^(t)E_ }dt 


(C.l) 


*'{(0  = C^e*!5  + C2afe  “lC  + J sinhtc^  U - t)  lAi^  (t)E+}dt 


£ 
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Similarly 


i£,(n)  = C5a2ea2r|  “ C6a2e  ^ + C~>  /cosh^a?(ri  " t)  }Ai{c,  (t)E+}dt 


7/c 


+ CgJcosh{a2(n  - t)  iAiit^ (t)E~ )dt 


:7a2  /. 


(C.4) 


i^Cn)  = Cj.a2ea2n  + C^a^e  a2n  + C7a2  I si-nll^cx 2^n  " t) }Ai{£2 (t)E+)dt 


C8a2  lsi' 

r.  * 


+ C^AK^C^E"*')  + Cg»2  / sinh{a2(n  - t)  } Ai { C2 (t)E~  }dt 


+ CgAi{;2(ir>)E-} 


(C.5) 


:7a2  Jc 


"^2  (h  ) = C^a|ea^n  - Cga|e  a2n  + j cosh{ (ft  - t)  }Ai{£2  (t)E+}dt 


+ C^Ai’  {^2  (n)E+}C2E+  + cga2  /cosh{a2(n  - t)  }Ai  { t,2  (t)E~  }dt 


+ CgAi' {c2(n)E-}^2E_ 


(C.6) 


where  it  follows  from  Eqs.  (4.3.4)  and  (4.3.20)  that 


A 
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APPENDIX  D 


DERIVATIVES  OF  GENERAL  SOLUTIONS  WITH  MASS  TRANSFER 


Differentiating  Eqs.  (5.4.15)  and  (5.5.16)  w.r.t.  4 and  Eqs. 
(5.4.17),  (5.5.22)  and  (5.5.27)  w.r.t.  n using  Leibniz's  rule;  and 

carrying  out  the  necessary  simplifications,  the  results  are 

4 

i|^(4)  = C1u1eal^  - C2a1e_ait’  + eRht//2  coshio^U  - t)  }Ai  Uj  (t) E+}dt 


* 


coshlat^  (4  - t) }Ai{q1(t)E~ }dt 


(D.l) 


■-a 

^(4)  = C1a^eal^  + C2a^e"al^  + c3a1  / eRh^2  sinMo^U  - t)  }Ai (c1  (t)E+)dt 

C 

+ C3eRh^^2  Ai(cl(OE+}  + c^a1l  eRht/2  sinhia^U  - t)  lAi^  (t)E+}dt 


+ C4eRh?/2  Ai(41(4)E-: 


(D.2) 


ip'"(0  = C]a3eal<’  - C2aJe-0llC  + C3a^eRht/2  coshUjU  - t)  >Ai{ Cj_  (t)E+}dt 


+ C3eV/2 


R. 

Ai{;i(4)E+}  + Ai'{;1(i;)E+}^E+ 


A 


+ C4l  / e h cosh{a1 (£  - t) }Ai{ci(t)E"}dt 


+ C4eRh^/2^  Al{Cl(OE-}  + Ai'U^OE-jqE-] 


(D.3) 


Similarly 

n 

*2^  = C5a2eCl2n  ” C6a2e  ^ + c7J^6t/2  sinh{u2(n  - t)  }Ai{  (t)E+}dt 

n* 

n 

+ Cg  /eR<5^2  sinh{a2(n  - t)  }Ai{^2  (t)E“}dt 


^2  (r|)  = C5a^ea2n 


+ C^a^e  a2n  + Cyd,,^* < 


(D.4) 


eR<$t^2  sinMo^Cn  - t)  }Ai(c2  (t)E+}dt 


/n 

eR6^/2  sinh{a2(ri  - t)  J Ai { C2  (t)E“  }dt 


+ C8eR6n/2  Ai{c2(n)E-} 


(D.5) 


'^2 Cn ) = Csa|ea2n  - Cfia^e  a2n  + cosh{a0(n  - t)  }Ai(c0(t)E+}dt 


:7eR6n/2  ~ Ai{;2(n)E+}  + Ai'  {^(rOE+J^E+j  + 
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+ C, 


rtf' 


|5t^2  cosh{a2(n  - t)  }Ai{;2(t)E"  }dt 


+ c8eR«n/2|^-2-  Aik2(n)E"}  + Ai'{;2(n)E-)e2E-] 

^(O  = eRhPrlC/2  ^ | C9Ai’{zi(OE+}E+  + C^Ai*  Z;L (^)E~  }E~  [* 

2itPr  R f|z’  ( I' 

+ 7,JClj;(0  + SJ2(°  + C3j3(°  + C4J4(° 

(o^r^Gp'  1 }- 


(D.6) 


+ eRhPrl^2 


V* 


~ I CQ Ai { z L ( ^ ) E+  j + C10Ai{Zl(C)E- 


+ 


2irPr  R T z I 

C1J1(^)  + c2J2(C)  + C3J3(C)  + C4J4(C) 

<VrxVi>1/3' 


(D.  7) 


where 


■/ 


RyPr, 

= / e(±al  2~~)l  G^tjOdt 


(D.8) 


■ lt-h 


J3,4U)  " /e'Rl>Prlt/2  j • 


^R^i/2  gjjjhfo,  _ t)  }Aik^(T)E±}dxdt 

(D.9) 


with 


G:(t;0  = Ai{z1(t)E+}Ai'{z1(C)E_}E~ 
-Ai{z1(t)E“}Ai'{z1(OE+}E+ 


(D . 10) 


z|  = -i(a1Pr1R1u1) 


Similarly 


§2  = eR<Sn/2  zj  C1;.Ai'{z2(n)E+}E+  + C^Ai ' { z2  (n)  E-  }E_ 


2^K2T  z l < . . 

+ jc^Cn)  + C6^6^n^  + C7^7^n^  + C8'’Vn'> 

x 1 / 31  ' j 

(a2R2u2 

+ eR6n/2^  C11Ai{z2(n)E+}  + C12Ai{z2(n)E“} 


2ttR2^2 


(a2R2u2 


fj/3f:5J5(Tl)  + C6J6(n)  + C7J7(n)  + C8J8(n)} 


^5,6(n)  =/*e(±a2  " 2 )C  G2(t;n)dt 


/>R^2  GAhrofW''2 


sinh{a„(t  - t)}A1(c„(t)E  }dxdt 


with 


G^( t;£)  * Ai{z2(t)E+}Ai' {z9 (n)E")E- 


-Ai{z2(t)Er}Ai'{z9(C)E+}E+ 


and 


'Z2  = ” 1 (a2R2^2^2 


Finally, 


x(n)  - e 


R5n/2 


,xj  C13Ai'{z2(n)E+}E+  + C14Ai'{z2(n)E-}E~j> 


2ttR2xz2 

, d ’ n 1 

(o^R^) 


;c5j5(n)  + c6J6(n)  + c737(o)  + cgJg(n)> 


+ eV/2 


R 3 

|c13Ai{z2(n)E+}  + c14Ai{z2(n)E-}| 


2^^ 


( ^2  ^2  u2  ^ 


i/3 1 5 5 


C1.j1.(n)  + c6J6(n)  + c?J7(n)  + c8Jg(n), 


(D.15) 


(D.16) 


(D. 17) 


with  Eqs . (D.13)  - (D.15)  applying  in  this  case  also. 
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APPENDIX  E 


LIST  OF  INTEGRALS  IN  ZERO  MASS  TRANSFER  PROBLEM 


The  integrals  1^  through  1^  in  Eqs.  (4.4.1)  through  (4.4.17) 


are  defined  as  follcws: 


-1 

IL  2 = ^sinh{a^(-l  - t)  }Ai{ £ ^(OE1  }dt 
-i 

I^  ^ = ycosh{a^(-l  - t)  }Ai(^^(t)E±)dt 


)}Ai{?2(t)E±idt 


1^  6 = ^sinh{ci2(l  - t)  lAit^^E^dt 

7 - 

l-j  g = /cosh{a2(l  - t)  }Ai{c2  (t)E' 

»*  0 

[9  10  =ysinh(a1t)Ai{;1(t)E1}dt 

£* 

^ 0 

12  ~ j cosh(a^t)Ai{i;^(t)E±}dt 


,u  -/• 


X13  14  = / 3inh(a2t)Ai{52(t)E±}dt 


Il5  16  = / cosh(a2t)Ai{^2 (t)E±}dt 
n*  * >’ 

with  and  q*  such  that  C^(£*)  = 0 and  ^2(n*)  = 0. 
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APPENDIX  F 


LIST  OF  INTEGRALS  IN  MASS  TRANSFER  PROBLEM 


The  integrals  1^  through  I^g  in  Eqs.  (5.6.1)  through  (5. 

are  defined  as  follows: 

-1 

1^  2 = / e^hc^  sinh{a^  (-1  - t)  }Ai{  (t) E1  }dt 

c*-i 

Ig  ^ cosh{a^(-l  - t)  }Ai(C^(t)E±}dt 

5*1 

1^  g = J'eP'6t^  sinh{a2(l  - t)  }Ai{ ^ (t) E1  }dt 
n*i 

= feR& 

7, S Je 

* 0 

I9  10  =f^/2  sinh(ait)Ai{c1(t)E±}dt 
r* 

... 


5 1/2 


I,  „ = /e‘<5  cosh{a2(l  “ t)  }Ai{  (t)E_  }dt 

* n 

n o 


Rht/2 


11 


cosh(a^t)Ai{  (t)E~}dt 


£*0 


In 


I^g  ^ sinh(a2t)Ai{C2 (t)E±}dt 

*0 

n/*  ' 

1^^  = /eR<5t/2  cosh(a2t)Ai{c2(t)E±}dt 

n* 

integrals  1^  ^ £ and  n are  such  that  £^(£  ) = 0 and 


1 


3 

6.15) 


(n*)  = 0. 


V 
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I17,18  = J 


I19,20  J 


1»2<_1) 

4*  _ 

iy-1’  -f( 


-1  . 
e(±al  " 2 )t:  F1(t;-l)dt 


Vrl 


1 2 


)t 


G1(t;-l)dt 


I21,22  J5,6 


r1  r6  ~ 

= Jc  *(1)  = /e(±a2  ' T>*  F2(t;l)dt 


R 


‘23,24  * -/e(i“2  ' 2“H  G2<t;1)dt 


I25,26  ~ J 


I27,23  = J 


1.2<0)  =f 
i-2<0)  ‘ J' 


n*0 

= |e(±al  " 2 


V*i)E 


F1(t;0)dt 


V1!, 


e(±al  - 2"  *)t  G^(t;0)dt 


29 ,30  °5,6 


6 

= Jc  *<0)  = l'.(±a2  ' 2~)t  F2(t;0)dt 


31,32  5,6 


f 

r\*J- 

Je 


n*  0 R6 

= Jc  ^(0)  = [e(±a 2 “ T"H  G2(t;0)dt 


In  integrals  I 7_32,  (,*  and  n*  are  such  that  z^S  ) = 0 and  z2(n*) 

F , F , G1  and  G„  are  given  by  Eqs.  (5.5.19),  (5.5.25),  (D.10)  and 
12  1 2 

(D . 15)  respectively . 

And  finally, 


*33,34  = J3 


I35,36  = J3 


4(-l)  = ^_RhPrlt/2  F1(t;-l)d1(t)dt 
n*  -i 

4(-l)  = le'R hPrlt/2  G1(t;-l)d1(t)dt 
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1 


/ 


1 

-Rrt/2 


X37  38  = J7,8(1)  = e 6 F2(t;-l)d2(t)dt 


nyi 

j7>8(1)  = /e"V/2  G2(t;-l)d2(t)dt 


I39,40  = J 


42  * J3  4(0)  ‘ _/i'Vrlt/2  F1(t;0)d1(t)dt 

5V°  - 

/e  RhPrlt/2  G1(t;0)d1(t)dt 

7,8(0)  = e~K&t/2  F2(t;0)d2(t)dt 

n*  0 

j7>8(°)  =ye-v 


*43,44  " J3,4(0) 


I45,46  = J 


I47,48  " J 


where 


and 


t 

i-L  (t)  =/eV 


G2(t;0)d2(t)dt 


sinh{oi^(t  - t) }Al{41(T)E±}dx 


d2(t) 


P 


/2 


sinh{a2(t  - t) }Ai(42 (x)E±}di 


(F.l) 


(F.  2) 


* * 

For  integrals  ^g,  £ , n and  t are  defined  such  that 

= ^(t*)  = 0 and  ^( t*)  - z0(t*)  = 0.  In  the  integrals 
1^  through  I^g!  4-^,  42»  z]_  ^d  z2  are  t*le  transFormations  given  by 
Eqs.  (5.4.8),  (5.4.18),  (5.5.3)  and  (5.4.18)  respectively.  These 
notations  used  in  the  mass  transfer  problem  are  the  same  as  in 
Chapter  IV  in  order  to  facilitate  comparisons. 
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APPENDIX  G 


EVALUATION  OF  SINGLE  INTEGRALS 


Consider  typical  single  integrals  encountered  in  the  problems 
with  and  without  mass  transfer  (Appendices  E and  F) 

5 


■/ 


.K't/2 


I =/e  h sinh{ot.  (C  - t)  }Ai{  (t)Ei  }dt 


(G.l) 


S* 
and 


I = 


= J 3inh{a^ 


(5  ~ t))Ai{c1(t)E±}dt 


(G.2) 


Eq.  (G.l)  reduces  to  (G.2)  when  = 0 and  hence  it  is  sufficient  to 
concentrate  on  Eq.  (G.l). 

A new  variable  of  integration  t^  defined  by  the  following 
equation  is  introduced: 


lit  ) = -S  +-S*-+  S-.i*  t 

lT  2 2 1 


(G.  3) 


Then  the  integral  (G.l)  can  be  written  as 
1 

' :i)/2 


lmLT^ 


-JeRht(tJ 


sinh  [a^{£;  - t^)}]  Ai  [ ^1(t(t1)}E±]  dt^ 

(G.4) 


1 


A 
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This  integral  is  in  a form  suitable  for  Gauss- Legendre  integration. 
In  fact,  Eq.  (G.4)  can  be  expressed  as 


I = 


_ C ~ £ 


-jfCtpdti 


(G.5) 


-1 


where 


f (t^)  = sinh  [o^U  - tCt^)}]  Ai  t (t^)  }E+  ] (G.6) 


Therefore, 


I = 


I - S 


wif(yi) 


(G.7) 


i=l 


where  are  the  zeros  of  Legendre  polynomials  and  w^  are  the 
corresponding  weight  factors.  n is  the  number  of  zeros.  The  function 
f in  Eq.  (G.6)  was  calculated  at  the  zeros  y^  using  an  Airy  function 
routine  described  in  Re.  45.  A modified  version  of  the  IBM  SSP-routine 
DQG32  was  employed  to  carry  out  the  summation  (G.7).  Solution  for 
the  final  eigenvalue  was  obtained  using  n = 16,  32  and  96.  It  was 
observed  that  the  difference  in  the  values  of  c^  with  32  and  96  points 
was  negligibly  small  and  therefore  a 32-point  scheme  was  adopted 
throughout  the  computational  work. 
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APPENDIX  H 


EVALUATION  OF  DOUBLE  INTEGRALS 


Consider  a typical  double  integral  in  Appendix  F: 

J,  i 


= je~Rh?rl^/2  F1(t;C) l J sinhla^t  - t) 


}Ai{;i(T)E+}drdt 


(H.l) 


Substituting  for  F1  from  Eq.  (5.5.19)  the  above  integral 


becomes 


I = Ai{z1(c)E~}I1  - Ai{Zl(OE+}I2 


(H.2) 


where 


Lx  2 = je~* hPrlt/2  Ai{z1(t)E±}J'eRliT 


^ sinh(a^ (t  - t) }Ai{ci(i)E+}dxdt 


(H.  3) 


By  writing  the  sinh  function  in  terms  of  exponentials  it  may  be 


verified  that 


I1  = 2 (S1AA  " S2AA^ 


(H.4) 


12  ~ 2 (S1BA  " S2BA) 


where 


/>  Khi'ri  . r c Rj, 

e<±al  ~ 2 Ai{Zl(t)E+}  / e(±a  + 2~)t  Ai{;1(T)E+}d^dt 


(H.  5) 


1BA.2BA 


/>t  RhPrl,-  f V~ 

e^'al  2 ^ Ai{z^(t)E-}  /e^_0t  + 2 Ai{c^(T)E+}dTdt 


(H.6) 


It  is  noticed  that  integrals  (H.5)  and  (H.6)  are  of  the  form 


£ t 


= Ja(  t)  Jl 

t* 


B(t ) drdt 


(H.  7) 


Thus  the  original  integral  (H.l)  has  been  expressed  in  terms  of 
four  simpler  iterated  integrals  S1M>  S2M>  SlgA  and  S^.  It  is 
possible  to  simplify  (H.7)  further  by  rewriting  it  in  the  form 


where 


S = / A(t)C(t)dt  = / D(t)dt 


C(t)  = / B(r)dT 


(H.  8) 


The  integral  (H.8)  can  be  evaluated  by  employing  the  transformation 
(G.3),  i.e. 
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^(tl)  ? 2 C1  + 


* 4 + C* 


Hence 


S = $■ 

b 2 


1 

-Ji){Ut1)}dt1 


-1 


It  then  follows  that 


i 

- j Mi 


S = ^ ^A{t(t1)}C{t(t1)}dt1 


Thus 


t(tx) 

Ctt^)  } =J B(t )dT 
t* 


The  following  transformation  is  now  introduced  analogous  to  E 


t^)  - t*  t(tL)  + t* 

*(t1>  = 2 T1  + 2 


Therefore  (H.13)  assumes  the  form 


C{  t (tl) J 


t(t.)  - t *r 

~2 1 B{T(t1)}dr1 


-1 


(H. 10) 
(H.ll) 

(H.12) 

(H.13) 

. (H. 10) , 
(H.14) 

(H. 15) 


J 
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Finally,  Eq.  (H.7)  reads 


£ - £*  / t^tl^ 

2 J 2 

-1 


£*  . [\ 

A{t(t1)}  }dT1dtL 


(H. 16) 


Eq.  (H. 16)  is  akin  to  Eq.  (G.4)  for  single  integrals  and  is  in  a 
form  convenient  for  Gauss- Legendre  integration  in  two  dimensions. 
Following  the  procedure  in  Ref.  46,  Eq.  (H.16)  can  be  put  in  the  form 


1 1 


€ - 6! 


‘hi  g(t1,T1)dT1c 


-1  -1 


1 1 

f y*h( t1,T1)dT1dt1 


(H. 17) 


(H. 18) 


t(t-i)  " t* 

h(tl,Tl)  = Hz A{t(t1)}B{i(T1) 


(H. 19) 


Hence 


J h(t1,T1)dT1  = K(tL) 


1 


1 

~J c(t1)dt1 


w k(y.) 

I 


(H. 21) 


where  y^  are  the  zeros  of  Legendre  polynomials  and  are  the 
corresponding  weight  factors,  n^  is  the  number  of  zeros  in  the 
interval  (-1,  1)  along  the  t^  axis. 


1 

=^h(yi,T1)dx1 


w . ,h(y . . . ) 
. lj  'l  'ij' 


(H . 22) 


where  y^  are  the  zeros  of  Legendre  polynomials  and  w_  are  the 
corresponding  weight  factors,  n ^ is  the  number  of  zeros  in  the 
interval  (-1,  1)  along  the  axis. 

The  final  form  of  the  integral  (H.16)  is 


e - v 
2 


nl  n2 


2>iZv(w 


(H. 23) 


i=l  j=l 


with  h defined  by  Eq.  (H.19). 

In  the  present  work,  experience  with  single  integrals  suggested 
the  choice  n^  = n^  = 32.  The  procedure  described  above  was  used  to 


r 
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calculate  the  integrals  S1AA>  S2AA’  S1BA’  anc*  S2BA*  Tll^s  computation 
leads  to  the  evaluation  of  1 through  equations  (H.4)  and  (H.2).  The 
details  of  the  computation  of  the  integrand  are  the  same  as  in 
Appendix  G. 


187 


APPENDIX  1 

NEWTON-RAPHSON  ITERATION  FOR  THE  EIGENVALUE 

As  mentioned  in  Secs.  4.5  and  5.7  the  eigenvalue  is 
evaluated  using  the  equations, 

(A(Cl))  {C>  = {V(Cl)}  (1.1) 

and 

G[c1,{C(c1)}]=  0 (1.2) 

For  the  zero  mass  transfer  problem,  [A(c^)]  is  an  8 x 8 coefficient 
matrix  of  the  left  hand  sides  of  Eqs.  (4.4.1)  - (4.4.8)  and  V(c^)is 
a column  matrix  of  the  right  hand  sides.  In  the  case  of  the  mass  trans- 
fer problem,  CA(c^)]  is  a 14  x 14  coefficient  matrix  of  the  left  hand 
sides  of  Eqs.  (5.6.1)  - (5.6.7)  and  (5.6.9)  - (5.6.15),  and  again 
V(c^)  is  a column  matrix  cf  the  right  hand  sides.  The  characteristic 
function  G is  given  by  Eq.  (4.4.11)  for  the  zero  mass  transfer  case 
and  by  Eq.  (5.6.17)  for  the  mass  transfer  problem. 


If  c^  is  a guess  value,  then  the  first  order  correction  due  to 


Newton-Raphson  method  is 


G(cx) 

AC1  = ' Gt(c^) 


(1.3) 


Thus  it  is  necessary  to  compute  the  derivative  G' (c^) . Differentia- 
ting (1.2)  w.r.t.  c^. 


r 


G'  (c,) 


3g_  V^3G_^i 

3c,  3C,  3c, 


(1.4) 


where  N = 8 for  the  zero  mass  transfer  and  N = 14  for  the  mass  transfer 
case.  The  calculation  of  3G/3c,  and  3G/3C,  from  Eqs.  (4.4.11)  and 
(5.6.17)  is  straight-forward.  The  calculation  of  3C,/3c,,  however, 
is  somewhat  involved  and  it  is  outlined  here. 


Differentiating  Eq,  (1.1)  w.r.t.  c,. 


[ A(cl)] 


3 { C } . 3[A(cl)] 
3Ci  " 3c, 


(C)  = 


3{V(Cl) } 


m = Zi  = rAr 

3c,  3c,  L -1 


(1.5) 


As  before,  {3V/3c,}  can  be  easily  obtained.  The  computation  of 
3[a]/3c,,  hcwever,  is  a very  complicated  task  because  it  invo 'ves 
differentiating  the  various  integrals  (in  Appendices  E and  F)  w.r.t. 
c,  using  Leibniz’s  rule.  This  leads  to  tedious  algebra  and  there- 
fore the  details  are  omitted.  It  is  sufficient  to  say  that  the 
derivatives  of  the  above-mentioned  integrals  can  be  related  to  the 
integrals  themselves  through  integration  by  parts.  Once  3[a]/3c, 
is  known  it  is  a simple  matter  to  compjte  G' (c,)  and  hence  Ac,. 


FIG.  I a.  STEADY-STATE  CONFIGURATION  WITHOUT  MASS 
TRANSFER 


FIG.  2 CONTINUITY  OF  TANGENTIAL  VELOCITIES  AT  INTERFACE 


FIG.  3 a.  BALANCE  OF  NORMAL  STRESSES  (zero  mass  transfer 
case) 


FIG.  3 b.  BALANCE  OF  NORMAL  STRESSES  (mass  transfer  case) 
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FIG.  4 FUNCTION  H(x)  VS  x 
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FIG.  7b.  PHASE  VELOCITY  CURVE  FOR  EIGENVALUE  C|2 


FIG.  9b.  PHASE  VELOCITY  CURVE  FOR  EIGENVALUE  C,4; 
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FIG.  12a.  COMPARISON  OF  AMPLIFICATION  CURVES  IN- 
CLUDING AND  EXCLUDING  INSTABILITIES  IN 
GAS  MOTION 
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FIG.  12b.  COMPARISON  OF  PHASE  VELOCITY  CURVES  IN- 
CLUDING AND  EXCLUDING  INSTABILITIES  IN 
GAS  MOTION 


mC2r/u'lO) 


m=  0 ) FOR  CONDITIONS 
0 | IN  TABLE  H 


